1 | ! **************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the subroutines: partem_s
|
---|
4 | !
|
---|
5 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
6 | ! Shura Hayryan, Chin-Ku
|
---|
7 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
8 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
9 | !
|
---|
10 | !
|
---|
11 | ! **************************************************************
|
---|
12 |
|
---|
13 | subroutine partem_s(num_rep, nequi, nswp, nmes, nsave, newsta,
|
---|
14 | & switch)
|
---|
15 | !
|
---|
16 | ! PURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM
|
---|
17 | ! ON SINGLE-PROCESSOR MACHINE.
|
---|
18 | !
|
---|
19 | ! CALLS: addang,energy,metropolis,rgyr,setvar,
|
---|
20 | !
|
---|
21 | include 'INCL.H'
|
---|
22 |
|
---|
23 | double precision gasc, dv, grnd, vr, addang, energ, energy, acz
|
---|
24 | double precision ee, rgy, temp0, delta, ra
|
---|
25 |
|
---|
26 | integer i, j, nml, k, nstart, iv, nsw, iswitch, k1, num_rep1, nu
|
---|
27 | integer in, jn
|
---|
28 |
|
---|
29 | external can_weight
|
---|
30 | ! TODO Store global coordinates in pgbpr
|
---|
31 | logical :: newsta
|
---|
32 | integer :: num_rep, nequi, nswp, nmes, nsave, switch
|
---|
33 | real*8 :: ttemp
|
---|
34 | character*80 :: filebase, fileNameMP
|
---|
35 |
|
---|
36 | parameter(gasc=0.00198773d0)
|
---|
37 | ! newsta: .true. for re-start; .false. for random start configurations
|
---|
38 | ! no: Number of replicas
|
---|
39 | ! nvrmax: Maximal number of dihdral angles
|
---|
40 | ! nequi: Number of sweeps for equilibrisation
|
---|
41 | ! nswp: Number of sweeps for simulation
|
---|
42 | ! nmes: Number of sweeps between measurements
|
---|
43 | ! gasc: Gas constant
|
---|
44 | !
|
---|
45 | real*8, allocatable :: coor_G(:, :),pbe(:),
|
---|
46 | & eol(:),acc(:),temp(:), pgbpr(:, :, :)
|
---|
47 | integer, allocatable :: ipoi(:)
|
---|
48 | ! temp: temperatures for each replica
|
---|
49 | ! coor_G: dihedral angles for ALL replicas
|
---|
50 | ! pbe: inverse temperatures for each replica
|
---|
51 | ! ipoi: Points to replica ipoi(k) which is currently
|
---|
52 | ! at inverse temperature pbe(k)
|
---|
53 | ! eol: energy of each replica
|
---|
54 | ! acc: accepatance rate of each replica
|
---|
55 | !
|
---|
56 | ! beta: inverse temperature of single replica
|
---|
57 | !
|
---|
58 | !
|
---|
59 | ! Allocate the arrays
|
---|
60 | allocate(coor_G(nvr, num_rep))
|
---|
61 | allocate(pgbpr(6, ntlml, num_rep))
|
---|
62 | allocate(temp(num_rep), pbe(num_rep), eol(num_rep))
|
---|
63 | allocate(acc(num_rep), ipoi(num_rep))
|
---|
64 |
|
---|
65 | if (.not.(allocated(coor_G).and.allocated(temp)
|
---|
66 | & .and. allocated(pbe).and. allocated(eol)
|
---|
67 | & .and. allocated(acc).and. allocated(ipoi))) then
|
---|
68 | stop "Unable to allocate memory in partem_s. Exiting."
|
---|
69 | end if
|
---|
70 | ! Initialize arrays
|
---|
71 | do i = 1, num_rep
|
---|
72 | do j = 1, nvr
|
---|
73 | coor_G(j, i) = 0.0
|
---|
74 | end do
|
---|
75 | do nml = 1, ntlml
|
---|
76 | do j = 1, 6
|
---|
77 | pgbpr(j, nml, i) = 0.0
|
---|
78 | end do
|
---|
79 | end do
|
---|
80 | temp(i) = 0.0
|
---|
81 | pbe(i) = 0.0
|
---|
82 | eol(i) = 0.0
|
---|
83 | acc(i) = 0.0
|
---|
84 | ipoi(i) = 0
|
---|
85 | end do
|
---|
86 |
|
---|
87 |
|
---|
88 | ! READ IN TEMPERATURES
|
---|
89 | open(11, file='temperatures', status='old')
|
---|
90 | do i=1,num_rep
|
---|
91 | read(11,*) j,ttemp
|
---|
92 | temp(j) = ttemp
|
---|
93 | pbe(j) = 1.0d0/(ttemp * gasc)
|
---|
94 | ipoi(j) = j
|
---|
95 | end do
|
---|
96 | close(11)
|
---|
97 |
|
---|
98 | open(13,file='time.d',status='unknown')
|
---|
99 |
|
---|
100 | if(.not.newsta) then
|
---|
101 | ! READ START Values
|
---|
102 | open(15,file='par_R.in',form='unformatted')
|
---|
103 | do k = 1, num_rep
|
---|
104 | read(15) nstart
|
---|
105 | read(15) eol(k)
|
---|
106 | lunvar = 14
|
---|
107 | filebase = 'conf_0000.var'
|
---|
108 | varfil = fileNameMP(filebase, 6, 9, k)
|
---|
109 | call redvar()
|
---|
110 | end do
|
---|
111 | close(15)
|
---|
112 | else
|
---|
113 | nstart=1
|
---|
114 | do i=1,num_rep
|
---|
115 | if (switch.ne.0) then
|
---|
116 | do j=1,nvr
|
---|
117 | iv = idvr(j)
|
---|
118 | if (switch.gt.0) then
|
---|
119 | dv=axvr(iv)*(grnd()-0.5)
|
---|
120 | vr=addang(pi,dv)
|
---|
121 | else
|
---|
122 | vr = pi
|
---|
123 | end if
|
---|
124 | coor_G(j,i)=vr
|
---|
125 | end do
|
---|
126 | end if
|
---|
127 | do nml = 1, ntlml
|
---|
128 | do j = 1, 6
|
---|
129 | pgbpr(j, nml, i) = gbpr(j, nml)
|
---|
130 | end do
|
---|
131 | end do
|
---|
132 | end do
|
---|
133 |
|
---|
134 | ! Equilibrization for each replica (No replica exchange move)
|
---|
135 | do k=1,num_rep
|
---|
136 | beta=pbe(k)
|
---|
137 | do i=1,nvr
|
---|
138 | iv =idvr(i)
|
---|
139 | vlvr(iv) = coor_G(i,k)
|
---|
140 | end do
|
---|
141 | do nml = 1, ntlml
|
---|
142 | do j = 1, 6
|
---|
143 | gbpr(j, nml) = pgbpr(j, nml, k)
|
---|
144 | end do
|
---|
145 | end do
|
---|
146 | energ = energy()
|
---|
147 | do nsw=1,nequi
|
---|
148 | CALL METROPOLIS(energ,acz,can_weight)
|
---|
149 | end do
|
---|
150 | write (logString, *)
|
---|
151 | & 'Start energy after equilibration for replica:',k, energ
|
---|
152 | do i=1,nvr
|
---|
153 | iv = idvr(i)
|
---|
154 | coor_G(i,k) = vlvr(iv)
|
---|
155 | end do
|
---|
156 | do nml = 1, ntlml
|
---|
157 | do j = 1, 6
|
---|
158 | pgbpr(j, nml, k) = gbpr(j, nml)
|
---|
159 | end do
|
---|
160 | end do
|
---|
161 | eol(k) = energ
|
---|
162 | end do
|
---|
163 | end if
|
---|
164 |
|
---|
165 | ! Now begins the simulation with Multiple Markov Chains
|
---|
166 | iswitch = 1
|
---|
167 | do nsw=nstart,nswp
|
---|
168 | !-------------------Ordinary MC sweep for each replica
|
---|
169 | do k1=1,num_rep
|
---|
170 | k = ipoi(k1)
|
---|
171 | beta = pbe(k1)
|
---|
172 | energ = eol(k)
|
---|
173 | acz= acc(k)
|
---|
174 | do i=1,nvr
|
---|
175 | iv = idvr(i)
|
---|
176 | vlvr(iv) = coor_G(i,k)
|
---|
177 | end do
|
---|
178 | do i=1,ntlml
|
---|
179 | do j = 1, 6
|
---|
180 | gbpr(j, nml) = pgbpr(j, i, k)
|
---|
181 | end do
|
---|
182 | call setvar(i,vlvr)
|
---|
183 | end do
|
---|
184 | CALL METROPOLIS(energ,acz,can_weight)
|
---|
185 | !
|
---|
186 | if(mod(nsw,nmes).eq.0) then
|
---|
187 | ! Measure and store here all quantities you want to analyse later
|
---|
188 | ! ee: end-to-end distance
|
---|
189 | ! rgy: radius of gyration
|
---|
190 | call rgyr(0,ee,rgy)
|
---|
191 | temp0=1.0d0/beta/gasc
|
---|
192 | write(13,'(i8,i3,4f12.3)') nsw,k1,temp0,energ,ee,rgy
|
---|
193 | end if
|
---|
194 | !
|
---|
195 | do i=1,nvr
|
---|
196 | iv = idvr(i)
|
---|
197 | coor_G(i,k) = vlvr(iv)
|
---|
198 | end do
|
---|
199 | do nml = 1, ntlml
|
---|
200 | do j = 1, 6
|
---|
201 | pgbpr(j, nml, k) = gbpr(j, nml)
|
---|
202 | end do
|
---|
203 | end do
|
---|
204 |
|
---|
205 | eol(k) = energ
|
---|
206 | acc(k) = acz
|
---|
207 | acz = 0.0d0
|
---|
208 | end do
|
---|
209 | !------------------Exchange of Replicas
|
---|
210 | if (.not.num_rep.gt.1) cycle ! Skip exchange if we only have a single replica.
|
---|
211 | if(iswitch.eq.1) then
|
---|
212 | num_rep1 = num_rep-1
|
---|
213 | nu = 1
|
---|
214 | iswitch = 2
|
---|
215 | else
|
---|
216 | num_rep1 = num_rep
|
---|
217 | nu = 2
|
---|
218 | iswitch = 1
|
---|
219 | end if
|
---|
220 | do i=nu,num_rep1,2 ! labels (inverse) temperatures
|
---|
221 | j=i+1
|
---|
222 | if(i.eq.num_rep) j=1
|
---|
223 | in=ipoi(i)
|
---|
224 | jn=ipoi(j)
|
---|
225 | delta=-pbe(i)*eol(jn)-pbe(j)*eol(in)
|
---|
226 | & +pbe(i)*eol(in)+pbe(j)*eol(jn)
|
---|
227 | ra = grnd()
|
---|
228 | if(ra.le.exp(delta)) then ! Metropolis to switch temperatures
|
---|
229 | ipoi(i) = jn ! between replicas
|
---|
230 | ipoi(j) = in
|
---|
231 | end if
|
---|
232 | end do
|
---|
233 | !-----------End exchange of Markov chains
|
---|
234 | !
|
---|
235 | end do
|
---|
236 | !
|
---|
237 | ! Write down final conformations for re-starts
|
---|
238 | open(15,file='par_R.in',form='unformatted')
|
---|
239 | do k=1,num_rep
|
---|
240 | write(15) nswp
|
---|
241 | write(15) eol
|
---|
242 | do i=1,nvr
|
---|
243 | iv = idvr(i)
|
---|
244 | vlvr(iv) = coor_G(i,k)
|
---|
245 | end do
|
---|
246 | do i=1,ntlml
|
---|
247 | do j = 1, 6
|
---|
248 | gbpr(j, nml) = pgbpr(j, i, k)
|
---|
249 | end do
|
---|
250 | call setvar(i,vlvr)
|
---|
251 | end do
|
---|
252 | filebase = "conf_0000.var"
|
---|
253 | call outvar(0, fileNameMP(filebase, 6, 9, k))
|
---|
254 | end do
|
---|
255 | close(15)
|
---|
256 | close(13)
|
---|
257 | return
|
---|
258 | end
|
---|
259 |
|
---|