-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain.py
151 lines (126 loc) · 7.97 KB
/
main.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
import numpy as np
import os
import time
from modules import initialize
from modules import dipole
from modules import computestatdc
from modules import correlationfunction
# ----------------------------------------------------------------------------------------------------------------------
# INITIALIZATION
start = time.time()
# INDICATES WHERE ARE THE DATA
root = '../dati10ns/'
filename = 'dump1.1fs.lammpstrj'
# GETS THE NUMBER OF PARTICLES IN THE SIMULATION
Npart = initialize.getNpart(filename, root)
# IF NOT ALREADY DONE SAVE THE DATA IN A BINARY FORMAT SO THAT IN THE FUTURE READS JUST THE BINARY
dati = initialize.getdatafromfile(filename, root, Npart) # JUST DO ONCE!!!! READ DATA FROM THE BINARY INSTEAD
# GETS THE DIMENSIONS OF THE SIMULATION BOX
Lato, Lmin = initialize.getBoxboundary(filename, root)
# GETS THE NUMBER OF SNAPSHOT OF THE SIMULATION. RESHAPE THE ARRAY OF THE DATA SO THAT WE HAVE A MATRIX WITH THE
# COORDINATES AND THE CHARGES FOR EACH SNAPSHOT
nsnapshot = initialize.getNsnap(dati, Npart)
print('The data are stored in '+root+filename)
print('The system has {}'.format(Npart)+' atoms in a box of side {:10.5f}'.format(Lato)+' Angstom')
print('In the calculation we are using {}'.format(nsnapshot)+' snapshots')
print('Initialization done in {:10.5f}'. format(time.time()-start)+'s')
# ----------------------------------------------------------------------------------------------------------------------
# CALCULATION OF THE DIPOLES
start1 = time.time()
# COMPUTES THE MATRIX OF THE MOLECULAR DIPOLES, THE CENTER OF MASS OF THE MOLECULE, THE ATOMIC CHARGES AND
# THE POSITION OF THE CHARGES (IN TIP4P/2005 THE OXY CHARGE IS IN A DIFFERENT POSITION THAN THE OXY ITSELF)
# POSO SETS THE DISTANCE BETWEEN THE OXY ATOM AND THE OXY CHARGE, IN THE TIP4P MODEL THE POSITIONS ARE DIFFERENT
posox = float(input('Set the OM distance for the calculation of the dipoles.>\n'))
dipmol, cdmol, chat, pos = dipole.computedipolestatdc(Npart, Lato, Lmin, nsnapshot, dati, posox)
print("Molecular dipoles, molecular positions, charges and charge positions for the trajectory computed in {:10.5f}".format(time.time()-start1)+'s')
# ----------------------------------------------------------------------------------------------------------------------
# CALCULATION OF THE STATIC DIELECTRIC CONSTANT
start2 = time.time()
nk = 120
# COMPUTES THE STATIC DIELECTRIC CONSTANT FOR NK VALUES OF THE G VECTOR IN THE (1,0,0) DIRECTION:
# 2\PI/LATO*(J,0,0), J=1,..NK
e0pol, e0ch = computestatdc.computestatdc(nk, dipmol, cdmol, chat, pos, Lato, 1)
e0pol, e0ch = computestatdc.computestatdc(nk, dipmol, cdmol, chat, pos, Lato, nsnapshot)
dipmol = []
cdmol = []
chat = []
pos = []
print('Static dielectric constant for {}'.format(nk)+' values of k computed in {:10.5f}'.format(time.time()-start2)+'s')
# ----------------------------------------------------------------------------------------------------------------------
# CALCULATION OF THE THERMOPOLARIZATION COEFFICIENT
start2 = time.time()
# COMPUTES THE THERMOPOLARIZATION COEFFICIENT FOR NK VALUES OF THE G VECTOR IN THE (1,0,0) DIRECTION:
# 2\PI/LATO*(J,0,0), J=1,..NK
chat, pos, enat, em, posatomic = dipole.computedipoletp(Npart, Lato, Lmin, nsnapshot, dati, posox)
tpc, sttpc = computestatdc.thermopolcoeff(nk, chat, enat, em, pos, posatomic, Lato, nsnapshot)
chat = []
pos = []
enat = []
em = []
posatomic = []
dipmol, cdmol, endip = dipole.computedipoletpdip(Npart, Lato, Lmin, nsnapshot, dati, posox)
tpcdip = computestatdc.thermopoldipcoeff(nk, dipmol, endip, cdmol, Lato, nsnapshot)
dipmol = []
cdmol = []
endip = []
print('Thermopolarization coefficient for {}'.format(nk)+' values of k computed in {:10.5f}'.format(time.time()-start2)+'s')
# ----------------------------------------------------------------------------------------------------------------------
# PRINT THE DIELECTRIC CONSTANT COMPUTED VIA THE POLARIZATION AND VIA THE CHARGES AS A FUNCTION OF G AND SAVE IN A FILE
file = '{}'.format(Npart)+'{}'.format(nk)+'dielconst.dat'
f = open(file, 'w+')
print("k\t"+'e0pol_xx\t'+'e0pol_yy\t'+'e0pol_zz\t'+'e0ch\n')
np.set_printoptions(precision=3)
f.write('#k (in units of 2pi/L(1,0,0))\t'+'e0pol_xx\t'+'e0pol_yy\t'+'e0pol_zz\t'+'e0ch\n')
for j in range(nk):
print('{}\t'.format(j)+'{:10.3f}\t'.format(e0pol[j][0])+'{:10.3f}\t'.format(e0pol[j][1])+'{:10.3f}\t'.format(e0pol[j][2])+'{:10.3f}'.format(e0ch[j]))
f.write('{}\t'.format(j)+'{:10.3f}\t'.format(e0pol[j][0])+'{:10.3f}\t'.format(e0pol[j][1])+'{:10.3f}\t'.format(e0pol[j][2])+'{:10.3f}\n'.format(e0ch[j]))
print('The static dielectric constants are saved in '+root+file)
print('The static dielectric constant for {}'.format(nk)+' values of k computed in : {:10.5f}'.format(time.time()-start2)+'s')
f.close()
# ----------------------------------------------------------------------------------------------------------------------
# PRINT THE DIELECTRIC CONSTANT COMPUTED VIA THE POLARIZATION AND VIA THE CHARGES AS A FUNCTION OF G AND SAVE IN A FILE
file = '{}'.format(Npart)+'{}'.format(nk)+'tpc.dat'
f = open(file, 'w+')
print("k\t"+'tpcd_xx\t'+'tpcd_yy\t'+'tpcd_zz\t'+'tpcch\t'+'std\n')
np.set_printoptions(precision=3)
f.write('#k (in units of 2pi/L(1,0,0))\t'+'tpcd_xx\t'+'tpcd_yy\t'+'tpcd_zz\t'+'tpcch\t'+'std\n')
for j in range(nk):
print('{}\t'.format(j) + '{:10.3f}\t'.format(tpcdip[j][0]) + '{:10.3f}\t'.format(tpcdip[j][1]) +\
'{:10.3f}\t'.format(tpcdip[j][2]) + '{:.2e}\t'.format(tpc[j]) + '{:.2e}'.format(sttpc[j]))
f.write('{}\t'.format(j) + '{:10.3f}\t'.format(np.real(tpcdip[j][0])) + '{:10.3f}\t'.format(np.real(tpcdip[j][1])) +\
'{:10.3f}\t'.format(np.real(tpcdip[j][2])) + '{:10.3f}\t'.format(np.real(tpc[j]))+ '{:10.5f}\n'.format(np.real(sttpc[j])))
print('The thermopolarization coefficient are saved in '+root+file)
print('The thermopolarization coefficient for {}'.format(nk)+' values of k computed in : {:10.5f}'.format(time.time()-start2)+'s')
f.close()
# ----------------------------------------------------------------------------------------------------------------------
# COMPUTES THE DIPOLES PAIR CORRELATION FUNCTION AT GMIN IN THE (G,0,0) DIRECTION AND SAVE IT IN A FILE
c = bool(int(input("Do you want to compute the dipole pair correlation function:True(1)/False(0)>\n")))
if c:
G = 0
start2 = time.time()
file = '{}'.format(Npart)+'{}'.format(G)+'dippcfwstd.dat'
print('The dipole pair correlation function for G=({}, 0, 0)'.format(G)+' is saved in '+root+file)
rdipmol, rcdmol, tdipmol, tcdmol = computestatdc.reshape(cdmol, dipmol)
nkk = 1
print('r\t'+'c_m_x\t'+'c_m_y\t'+'c_m_z\t'+'std_c_m_x\t'+'std_c_m_y\t'+'std_c_m_z\t')
gk, stdgk = computestatdc.dip_paircf(G, nkk, rdipmol, rcdmol, tdipmol, tcdmol, Lato, nsnapshot)
print(time.time()-start2)
start2=time.time()
nkk = 100
print('r\t'+'c_m_x\t'+'c_m_y\t'+'c_m_z\t'+'std_c_m_x\t'+'std_c_m_y\t'+'std_c_m_z\t')
gk, stdgk = computestatdc.dip_paircf(G, nkk, rdipmol, rcdmol, tdipmol, tcdmol, Lato, nsnapshot)
f = open(file, 'w+')
for i in range(nk):
f.write('{:10.5f}\t'.format((i*(Lato - 2)/2/nkk + 2)/0.529) + '{:10.5f}\t'.format(gk[i][0]) + '{:10.5f}\t'.format(gk[i][1]) + '{:10.5f}\t'.format(gk[i][2])+\
'{:10.5f}\t'.format(stdgk[i][0]) + '{:10.5f}\t'.format(stdgk[i][1]) + '{:10.5f}\n'.format(stdgk[i][2]))
print('The dipole pair correlation function for G=({}, 0, 0)'.format(G)+' computed in : {:10.5f}'.format(time.time()-start2)+'s')
print('Total elapsed time: {:10.5f}'.format(time.time()-start)+'s')
f.close()
# ----------------------------------------------------------------------------------------------------------------------
# COMPUTES THE DENSITY-DENSITY CORRELATION FUNCTION AS A FUNCTION OF K AND \OMEGA AND PRINTS IT IN A FILE
enat, posatomic = dipole.computedipolecorren(Npart, Lato, Lmin, nsnapshot, dati, posox)
nk = int(input('How many k points do you want in the density-density correlation function>'))
ftaf = correlationfunction.correlation(nk, nsnapshot, Lato, enat, posatomic)
for i in range(ftaf.shape[0]):
for j in range(ftaf.shape[1]):
pass