-
Notifications
You must be signed in to change notification settings - Fork 8
/
SMD_to_FE.f90
185 lines (144 loc) · 5.55 KB
/
SMD_to_FE.f90
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!@ @
!@ This program uses Jarzynski relationship to calculate free energy @
!@ profiles from work done in independent SMD simulations @
!@ @
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Program written by Jason Swails, jason.swails@gmail.com
! Last update 7/23/2008
!
!
! Description: This program uses the Jarzynski relationship to convert a collection of
! work data compiled from multiple steered molecular dynamics simulations into a free
! energy profile for the process. I have a script that complements this program,
! compiling a single file in the format read by this program from raw data files generated
! by AMBER's SMD module.
subroutine get_num_tokens(string, token_num)
implicit none
! Passed arguments
character(*), intent(in) :: string
integer, intent(out) :: token_num
! Local variables
integer :: string_loc ! our location in the string
integer :: iend ! last non-whitespace character location
string_loc = 1
iend = len_trim(string)
token_num = 0
do while (string_loc .le. iend)
if ( string(string_loc:string_loc) .le. ' ' ) then
string_loc = string_loc + 1
else
do while ( string(string_loc:string_loc) .gt. ' ' )
string_loc = string_loc + 1
end do
token_num = token_num + 1
end if
end do
end subroutine get_num_tokens
program SMD_to_FE
implicit none
real, dimension(:), allocatable :: Work_array
real, parameter :: KB = 1.98722E-3
real :: TEMP, distance, BETA, work, holder
integer :: num_trials
character (len=256) :: dataFile, outputFile
character (len=2048) :: line
character (len=3) :: response
character (len=1) :: input_open
integer :: i, c
! Work_array stores the values of work for each simulation: updated every new distance
! KB = Boltzmann's constant in kcal/mol*K (R, or gas constant)
! TEMP = system temperature, input by the user and checked for mistakes
! distance = the distance at which the Work values occur in the file
! num_trials = number of trials
! dataFile = input name for the file containing the data
! outputFile = name of file created to store output
! response = yes/no response to queries
! input_open = indicator so that error message is not displayed before error is made
! i, c = counters used in loops to keep track of indices
! BETA = stat-mech beta, 1/kT
! work = keeps track of exponential of work
! holder = holding variable that holds a variable's value while it is changed
!==================================================================!
! Get the temperature and impose checks to make sure it was entered correctly, i.e. Kelvins
! and no obvious typos made
!============================ GET TEMPERATURE =============================================
do
print*, "Enter the temperature of the system."
read (*,*) TEMP
if (TEMP < 0) then
print*, "Invalid Temperature. Enter system temp in K"
else if (TEMP < 100) then
write (*,*) "Is the temperature really ", TEMP, " K?"
read (*,*) response
if (response == "yes" .or. response == "Yes" .or. response == "YES" .or. response == "y" &
.or. response == "Y" .or. response == "YEs" .or. response == "yES") then
exit
end if ! response == ...
else if (TEMP > 1000) then
write (*,*) "Is the temperature really ", TEMP, " K?"
read (*,*) response
if (response == "yes" .or. response == "Yes" .or. response == "YES" .or. response == "y" &
.or. response == "Y" .or. response == "YEs" .or. response == "yES") then
exit
end if ! response == ...
else
exit
end if ! TEMP <, etc.
end do
!========================= END OF GET TEMPERATURE ========================================
! Ask for the input file, and check to make sure it exists. User can exit by typing
! quit here or in the output file prompt, which comes right after
input_open = "n"
10 if ( input_open == "y" ) then
print*, "There was an error opening the file, enter another one."
end if
print*,"Enter the name of the file with the data."
read(*,*) dataFile
if ( dataFile == "quit" ) then
goto 12
end if
input_open = "y"
open(7,file=dataFile,STATUS='OLD',ERR=10)
print*,"Enter the desired name of the output file."
read(*,*) outputFile
if( outputFile == "quit" ) then
goto 12
end if
!======================== END OF GET FILE NAMES ========================================
! Here, determine the number of trials that were done
i = 1
read(7, fmt='(a)') line
call get_num_tokens(line, num_trials)
rewind(7)
!======================== END OF TRIAL NUM DETERMINATION ===============================
! Set initial values
work = 0
BETA = 1 / (KB * TEMP)
allocate (Work_array(num_trials))
! Create the output file
open(8, file=outputFile)
! Implement Jarzynski's method
! exp(FE(x)/kT) = < exp(W(x)/kT) >
do
read(7,*,end=11)(Work_array(i), i=1,num_trials)
distance = Work_array(1)
do c=2, num_trials
work = exp(-BETA * Work_array(c)) + work
end do
holder = - log (work / (num_trials - 1) ) / BETA
write (8,*) distance, holder
work = 0
end do
! End the program, deallocate the array and close the files
11 continue
deallocate(Work_array)
close(7)
close(8)
print*, "Your free energy profile is complete! ", outputFile
12 continue
if( dataFile == "quit" .or. outputFile == "quit") then
print*, "Program terminated. Nothing computed."
end if
end program SMD_to_FE