[e40e335] | 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 | c Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
| 8 | c 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 |
|
---|