Scripts_for_NatComm/ion_time_series.inp at main · hzhekova/Scripts_for_NatComm · GitHub

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
* GENERATED BY CHARMM-GUI (http://www.charmm-gui.org)
* SETUP PERIODIC BOUNDARY CONDITION
*

!This CHARMM script evaluates the minimum distances between specified ions and protein atoms, or between specified ions at every step of a MD trajectory (the topology and trajectory files are in NAMD format) and have been used to generate the ion time series plots, presented in Fig.S5-S7 of the manuscript


bomblev -2
DIMENS CHSIZE 3000000 MAXRES 3000000

!Stream the necessary topology and parameter files
stream ../toppar.str


!Open the psf and dcd files from a MD simulation and the output file name_of_file_mind_coord.dat
open read unit 10 card name name_of_psf.psf
read psf  xplo unit 10 card

open unit 32 write form name name_of_file_mind_coord.dat
write title unit 32
*#Time  CO3-S1   Na-S1  Na1-CO3 CL-S1
*

open unit 21 read unform name name_of_dcd.dcd

!Read the DCD trajectory. Here the trajectory is loaded up to the 5200th frame
traj iread 21 nunits 1 begin 0 skip 1 stop 5200


!Define the protein atoms that will serve as a measure of distance from site S1
define s1 sele segid PROA .and. resid 846 .and. type CA end

!Define the anion, for which the time series are calculated
define co3 sele segid HETA .and. type C* end


!Start looping through the trajectory and for each frame evaluate the minimum distance from the CO3 ion to sites S1 (CO3-S1), the minimum distance from the closest Na ion to site S1, the minimum distance between the closest to site S1 Na ion and the CO3 ion (Na1-CO3), and the m,inimum distance between the closest CL and site S1.
set i = 1
label loop
traj read

coor mind sele co3 end sele s1 end
set co3s1 = ?mind

coor mind sele type SOD end sele s1 end
set s1dist = ?mind
set sod1 ?minda1  !Pick up Na+ #
coor mind sele bynu @sod1 end sele co3 end
set na1co3 ?mind  !Pick up Na+ #

coor mind sele type CLA end sele s1 end
set chlorid1 ?mind

#Convert the trajectory steps to ns (the conversion factor depends on the trajectory saving settings during the MD simulations)
calc time = @i * 0.24

coor stat
coor stat sele segid MEM* end
coor stat sele segid PRO* end


!Write  the information in the name_of_file_mind_coord.dat file
write title unit 32
*@time @co3s1 @s1dist @na1co3 @chlorid1


!Continue looping until the last trajectory step (in this case 5200th) is reached
incr i by 1
if i lt 5200 goto loop


stop