-
Notifications
You must be signed in to change notification settings - Fork 0
/
pdb.main_nongauss.py
42 lines (33 loc) · 1.68 KB
/
pdb.main_nongauss.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
from class_data import data
if __name__ == '__main__':
import time as timer
pdbfilename = '../nongauss_every100ps.pdb'
##--- output general settings ---=
import numpy as np
fn_prefix = '../AmC2_new_'
##--- timer start ---
start = timer.perf_counter()
##--- initialize a data object ---
d1 = data()
d1.save_mem=1
##--- read in pdb traj output ---
AN_gen_kw = [ range(11,11+400), 1, 1, 1, [], 15, "sel[:]", 'type']
# CT_gen_kw = [ range(1,1+10), 2, 40, 1, [30,1171], 15, "sel[ ((sel.index%30<=9 ) & (sel.index%30>=5)) | ( (sel.index%30>=15 ) & (sel.index%30<=24)) ]", 'type' ]
d1.read_all_pdb(pdbfilename, AN_gen = AN_gen_kw )#, CT_gen =CT_gen_kw )
#print(d1.allframes[0].L_AN)
print('last timestep read: ', d1.allframes[-1].time)
##--- unwrap anion coordinates --- # may not necessary for some cases
#d1.unwrapall_AN() # do not unwrap if read from pdb. pdb has unwrapped data already.
#print('unwrapped data')
##--- non gaussian parameter ---
time_col, nongauss_col = d1.nongauss_AN_avg(100, 800) # start_interval = 1 frame, 100ps/frame, maxattemp =500 by default
print('non gauss calculated')
##--- MSD
time1_col, msd_col = d1.msd_AN_avg(100, 800)
print('MSD calculated')
##--- timer stop ---
stop = timer.perf_counter()
print('time used in sec: ', stop-start)
##--- save non gaussian parameter
np.savetxt( fn_prefix+'nongauss.dat', np.transpose( [time_col, nongauss_col] ), fmt=['%d', '%f'], header='time/ps, nongauss' )
np.savetxt( fn_prefix+'MSD.dat', np.transpose( [time1_col, msd_col] ), fmt=['%d', '%f'], header='time/ps, msd' )