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