[e40e335] | 1 | ! **************************************************************
|
---|
| 2 | !
|
---|
| 3 | ! This file contains the subroutines: partem_s
|
---|
| 4 | !
|
---|
| 5 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
[32289cd] | 6 | ! Shura Hayryan, Chin-Ku
|
---|
[bd2278d] | 7 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
| 8 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
[e40e335] | 9 | !
|
---|
| 10 | !
|
---|
| 11 | ! **************************************************************
|
---|
[32289cd] | 12 |
|
---|
[e40e335] | 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'
|
---|
[32289cd] | 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 |
|
---|
[e40e335] | 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)
|
---|
[32289cd] | 53 | ! eol: energy of each replica
|
---|
[e40e335] | 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))
|
---|
[32289cd] | 64 |
|
---|
[e40e335] | 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
|
---|
[32289cd] | 86 |
|
---|
[e40e335] | 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')
|
---|
[32289cd] | 99 |
|
---|
[e40e335] | 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
|
---|
[32289cd] | 131 | end do
|
---|
[e40e335] | 132 | end do
|
---|
[32289cd] | 133 |
|
---|
[e40e335] | 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
|
---|
[38d77eb] | 150 | write (logString, *)
|
---|
| 151 | & 'Start energy after equilibration for replica:',k, energ
|
---|
[e40e335] | 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
|
---|
[32289cd] | 164 |
|
---|
[e40e335] | 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)
|
---|
[32289cd] | 185 | !
|
---|
[e40e335] | 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
|
---|
[32289cd] | 204 |
|
---|
[e40e335] | 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
|
---|
[32289cd] | 223 | in=ipoi(i)
|
---|
[e40e335] | 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 |
|
---|