Changeset 32289cd
- Timestamp:
- 11/19/09 11:29:41 (14 years ago)
- Branches:
- master
- Children:
- 38d77eb
- Parents:
- 6650a56
- Files:
-
- 35 edited
Legend:
- Unmodified
- Added
- Removed
-
INCL.H
r6650a56 r32289cd 1 ! implicit none 2 implicit integer*4 (i-n) 3 implicit real*8 (a-h,o-z) 1 implicit none 2 4 3 character*255 version 5 4 -
bgs.f
r6650a56 r32289cd 2 2 ! Trial version implementing the semi-local conformational update 3 3 ! BGS (Biased Gaussian Steps). This file presently contains the 4 ! functions initlund, bgsposs and bgs. 4 ! functions initlund, bgsposs and bgs. 5 5 ! 6 6 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, … … 9 9 10 10 11 ! Checks if it is possible to perform a BGS update starting at the 11 ! Checks if it is possible to perform a BGS update starting at the 12 12 ! variable indexed ipos. Calls: none. 13 13 logical function bgsposs(ips) 14 14 include 'INCL.H' 15 15 include 'incl_lund.h' 16 integer jv, ips, iaa, nursvr, nnonfx, i 17 16 18 logical ians 17 19 … … 20 22 ians=.true. 21 23 ! print *,'evaluating bgs possibility for ',ips,nmvr(jv) 22 if (nmvr(jv).ne.'phi') then 24 if (nmvr(jv).ne.'phi') then 23 25 ! print *,'bgs not possible because variable name is ',nmvr(jv) 24 26 ians=.false. 25 else if (iaa.gt.(irsml2(mlvr(jv))-3)) then 27 else if (iaa.gt.(irsml2(mlvr(jv))-3)) then 26 28 ! print *,'bgs impossible, residue too close to end' 27 29 ! print *,'iaa = ',iaa,' end = ',irsml2(mlvr(jv)) 28 30 ians=.false. 29 else 30 nnonfx=0 31 else 32 nnonfx=0 31 33 do i=iaa,iaa+3 32 if (iphi(i).gt.0) then 34 if (iphi(i).gt.0) then 33 35 if (.not.fxvr(iphi(i))) then 34 36 nnonfx=nnonfx+1 35 37 endif 36 38 endif 37 if (.not.fxvr(ipsi(i))) then 39 if (.not.fxvr(ipsi(i))) then 38 40 nnonfx=nnonfx+1 39 41 endif … … 52 54 53 55 ! Biased Gaussian Steps. Implements a semi-local conformational update 54 ! which modifies the protein backbone locally in a certain range of 55 ! amino acids. The 'down-stream' parts of the molecule outside the 56 ! which modifies the protein backbone locally in a certain range of 57 ! amino acids. The 'down-stream' parts of the molecule outside the 56 58 ! region of update get small rigid body translations and rotations. 57 59 ! 58 ! Use the update sparingly. It is rather local, and is not of great 60 ! Use the update sparingly. It is rather local, and is not of great 59 61 ! value if we need big changes in the conformation. It is recommended 60 ! that this update be used to refine a structure around a low energy 62 ! that this update be used to refine a structure around a low energy 61 63 ! candidate structure. Even at low energies, if you always 62 ! perform BGS, the chances of coming out of that minimum are small. 63 ! So, there is a probability bgsprob, which decides whether BGS or the 64 ! normal single angle update is used. 64 ! perform BGS, the chances of coming out of that minimum are small. 65 ! So, there is a probability bgsprob, which decides whether BGS or the 66 ! normal single angle update is used. 65 67 ! 66 68 ! Calls: energy, dummy (function provided as argument), addang, (rand) … … 70 72 include 'INCL.H' 71 73 include 'incl_lund.h' 74 double precision grnd, xiv, bv, ab, rv, dv, a, sum, p, r1, r2 75 double precision ppsi, wfw, addang, enw, energy, wbw, rd, delta 76 double precision dummy, eol1 77 78 integer ivar, i, nph, jv, ia, nursvr, icurraa, j, k, l 79 72 80 external dummy 73 81 dimension xiv(8,3),bv(8,3),rv(3,3),dv(3,8,3) … … 76 84 ! Initialize 77 85 ! print *,'using BGS on angle ',nmvr(idvr(ivar)) 78 if (bgsnvar.eq.0) then 86 if (bgsnvar.eq.0) then 79 87 bgs=0 80 88 goto 171 … … 91 99 do i=1,4 92 100 icurraa=ia+i-1 93 if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 101 if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 94 102 nph=nph+1 95 103 xiv(nph,1)=xat(iCa(icurraa)) … … 102 110 & +bv(nph,3)*bv(nph,3) 103 111 iph(nph)=iphi(icurraa) 104 endif 112 endif 105 113 if (.not.fxvr(ipsi(icurraa))) then 106 114 nph=nph+1 … … 136 144 enddo 137 145 enddo 138 do i=1,nph 146 do i=1,nph 139 147 do j=i,nph 140 148 A(i,j)=0 … … 157 165 sum=sum-A(i,k)*A(j,k) 158 166 enddo 159 if (i.eq.j) then 160 p(i)=sqrt(sum) 161 else 167 if (i.eq.j) then 168 p(i)=sqrt(sum) 169 else 162 170 A(j,i)=sum/p(i) 163 171 endif … … 177 185 dph(i)=0 178 186 enddo 179 ! Solve lower triangular matrix to get dphi proposals 187 ! Solve lower triangular matrix to get dphi proposals 180 188 do i=nph,1,-1 181 189 sum=ppsi(i) … … 190 198 do i=1,nph 191 199 sum=sum+ppsi(i)*ppsi(i) 192 enddo 200 enddo 193 201 wfw=exp(-sum) 194 202 do i=1,nph … … 208 216 do i=1,4 209 217 icurraa=ia+i-1 210 if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 218 if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 211 219 nph=nph+1 212 220 xiv(nph,1)=xat(iCa(icurraa)) … … 219 227 & +bv(nph,3)*bv(nph,3) 220 228 iph(nph)=iphi(icurraa) 221 endif 222 if (.not.fxvr(ipsi(icurraa))) then 229 endif 230 if (.not.fxvr(ipsi(icurraa))) then 223 231 nph=nph+1 224 232 xiv(nph,1)=xat(iC(icurraa)) … … 253 261 enddo 254 262 enddo 255 do i=1,nph 263 do i=1,nph 256 264 do j=i,nph 257 265 A(i,j)=0 … … 262 270 enddo 263 271 A(i,j)=bbgs*A(i,j) 264 if (i.eq.j) then 272 if (i.eq.j) then 265 273 A(i,j)=A(i,j)+1 266 274 endif … … 274 282 sum=sum-A(i,k)*A(j,k) 275 283 enddo 276 if (i.eq.j) then 277 p(i)=sqrt(sum) 278 else 284 if (i.eq.j) then 285 p(i)=sqrt(sum) 286 else 279 287 A(j,i)=sum/p(i) 280 288 endif … … 285 293 do j=i+1,nph 286 294 ppsi(i)=ppsi(i)+A(j,i)*dph(j) 287 enddo 295 enddo 288 296 enddo 289 297 sum=0 … … 291 299 sum=sum+ppsi(i)*ppsi(i) 292 300 enddo 293 wbw=exp(-sum) 301 wbw=exp(-sum) 294 302 do i=1,nph 295 303 wbw=wbw*p(i) … … 306 314 ! call outpdb(0,'after.pdb') 307 315 ! print *,'after outpdb for after.pdb' 308 ! do i=1,nph 316 ! do i=1,nph 309 317 ! print *,'BGS>',i,iph(i),vlvr(iph(i)),dph(i) 310 318 ! enddo 311 if (rd.ge.delta) then 319 if (rd.ge.delta) then 312 320 ! accept 313 321 eol1=enw 314 322 bgs=1 315 323 ! print *,'BGS move accepted' 316 else 324 else 317 325 ! reject 318 326 vlvr = ovr -
cnteny.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 20 20 include 'INCL.H' 21 21 22 parameter (coeycn=2.d0) ! min. cont. energy to display 22 double precision coeycn, coey, eyatcn, cqi, xi, yi, zi, xij, yij 23 double precision zij, rij2, rij4, rij6, rij, sr, ep, ey 24 25 integer ieyel, ntlvr, nml, iat1, nat, i, ifivr, i1s, io, iv, ia 26 integer it, ic, i2s, ims, i1, i2, ii, ity, ivw, j, jj, jty, i14 27 integer nbc, ir, nursat, idxat 28 29 parameter (coeycn=2.d0) ! min. cont. energy to display 23 30 24 31 dimension eyatcn(mxat),idxat(mxat) -
difang.f
r6650a56 r32289cd 21 21 ! ...................................................... 22 22 23 implicit real*8 (a-h,o-z)24 23 implicit none 24 double precision pi, pi2, d, a2, a1 25 25 parameter (pi=3.141592653589793d0, 26 26 & pi2=2.d0*pi) … … 47 47 ! ...................................................... 48 48 49 implicit real*8 (a-h,o-z) 49 implicit none 50 double precision pi, pi2, d, a1, a2 50 51 51 52 parameter (pi=3.141592653589793d0, -
energy.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 22 22 23 23 include 'INCL.H' 24 double precision esm, teysl, enyshe, enyflx, enylun, enyreg 25 double precision enysol, esolan, exvlun, eyabgn, eninteract 26 27 integer i 28 24 29 double precision teysm, teyel, teyvw, teyhb, teyvr 25 ! print *,'energy function with ientyp = ',ientyp 30 ! print *,'energy function with ientyp = ',ientyp 26 31 esm = 0.d0 27 32 teysm = 0.d0 … … 32 37 teysl = 0.d0 33 38 34 do i = 1,ntlml 39 do i = 1,ntlml 35 40 eysm=0 36 41 eyel=0 … … 44 49 if (ientyp.eq.0.or.ientyp.eq.3) then 45 50 esm=esm+enyshe(i) 46 else if (ientyp.eq.1) then 51 else if (ientyp.eq.1) then 47 52 esm=esm+enyflx(i) 48 53 else if (ientyp.eq.2) then … … 57 62 ! The Lund term stores the hydrophobicity energy in eysl 58 63 teysl = teysl + eysl 59 else 64 else 60 65 ! .. and the excluded volume term in eyvw, which is calculated once. 61 66 teyvw = teyvw + eyvw 62 67 endif 63 68 64 if (ireg.eq.1) eyrg=enyreg(i) 69 if (ireg.eq.1) eyrg=enyreg(i) 65 70 66 71 enddo 67 72 68 if (ientyp.ne.2) then 69 ! Don't touch eysl if using Lund potential, as enylun stores 73 if (ientyp.ne.2) then 74 ! Don't touch eysl if using Lund potential, as enylun stores 70 75 ! its hydrophobicity term there. 71 76 if (itysol.gt.0) then … … 79 84 eysl=0.d0 80 85 endif 81 else 86 else 82 87 ! Add excluded volume term and save it in eyvw 83 88 esm=esm+exvlun(0) … … 85 90 endif 86 91 87 ! The Abagyan entropic corrections depend on the area exposed to the 92 ! The Abagyan entropic corrections depend on the area exposed to the 88 93 ! solvent for each residue. So, this term has to be evaluated after the 89 94 ! solvent term. 90 95 eyab=0.0 91 96 if (ientyp.eq.3) then 92 do i = 1,ntlml 97 do i = 1,ntlml 93 98 eyab=eyab+eyabgn(i) 94 99 enddo … … 103 108 eyvr = teyvr 104 109 eysl = teysl 105 110 106 111 if (ientyp.ne.2) then 107 112 ! This is temporary. eninteract() does not yet know how to calculate … … 116 121 !c Calculates the internal energy for a single molecule. 117 122 ! All the partial energies are thus set to their values for molecule 118 ! nml. 123 ! nml. 119 124 ! 120 125 ! @param nml the ID of the molecule … … 125 130 126 131 !f2py intent(in) nml 127 132 128 133 include 'INCL.H' 134 double precision esm, enyshe, enyflx, enylun, enyreg, enysol 135 double precision esolan, exvlun, eyabgn 136 137 integer nml, i 138 129 139 esm = 0.d0 130 140 131 call setvar(nml,vlvr) 141 call setvar(nml,vlvr) 132 142 133 143 if (ientyp.eq.0.or.ientyp.eq.3) then … … 149 159 eysl=0.d0 150 160 endif 151 else 161 else 152 162 esm=esm+exvlun(nml) 153 163 endif 154 ! The Abagyan entropic corrections depend on the area exposed to the 164 ! The Abagyan entropic corrections depend on the area exposed to the 155 165 ! solvent for each residue. So, this term has to be evaluated after the 156 166 ! solvent term. -
eninteract.f
r6650a56 r32289cd 15 15 16 16 include 'INCL.H' 17 17 18 double precision cqi, xi, yi, zi, xj, yj, zj, xij, yij, zij, rij2 19 double precision rij4, rij6, rij, sr, ep 20 21 integer iml, jml, ires, jres, iat, ity, jat, jty 22 18 23 eysmi=0.0 19 24 eyeli=0.0 … … 42 47 yj=yat(jat) 43 48 zj=zat(jat) 44 49 45 50 xij=xat(jat)-xi 46 51 yij=yat(jat)-yi … … 69 74 & -chb(ity,jty)/(rij6*rij4) 70 75 endif 71 enddo 76 enddo 72 77 enddo 73 78 enddo … … 78 83 eysmi = eyeli + eyvwi + eyhbi 79 84 eninteract = eysmi 80 return 81 end 85 return 86 end -
enylun.f
r6650a56 r32289cd 1 1 ! ******************************************************************* 2 2 ! SMMP version of Anders Irback's force field, to be called the Lund 3 ! force field. This file contains the function enylun, which in turn 4 ! calls all the terms in the energy function. The terms Bias (ebias), 5 ! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the 6 ! Excluded volume (eexvol and eloexv) are also implemented in this 3 ! force field. This file contains the function enylun, which in turn 4 ! calls all the terms in the energy function. The terms Bias (ebias), 5 ! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the 6 ! Excluded volume (eexvol and eloexv) are also implemented in this 7 7 ! file. 8 8 ! … … 13 13 include 'INCL.H' 14 14 include 'incl_lund.h' 15 integer iat1, iat2 16 15 17 integer ires,frst,lst,root 16 18 do iat1=iCa(ires)+frst,iCa(ires)+lst … … 30 32 include 'INCL.H' 31 33 include 'incl_lund.h' 34 double precision eunit, csacc 35 36 integer i, j, iml, iat1, iat2, k, iat3, l, iat4, m, iat3p, irs 37 integer iatoff, iatmrg, ilp, lci 38 32 39 character mynm*4 33 40 logical prlvr … … 83 90 prlvr=.true. ! residue i is a proline variant 84 91 print *, 'proline variant ',mynm,i 85 else 92 else 86 93 prlvr=.false. 87 endif 88 89 if (mynm.eq.'ala') then 94 endif 95 96 if (mynm.eq.'ala') then 90 97 nhpat(i)=1 91 98 ihpat(i,1)=iCa(i)+2 … … 193 200 194 201 ! Initialization of the connections matrix matcon(i,j). The index 195 ! i runs from -mxconr to +mxconr, and j from 1 to mxat. 202 ! i runs from -mxconr to +mxconr, and j from 1 to mxat. 196 203 ! matcon(i2-i1,i1) = 0, if the distance between atoms i1 and i2 is fixed 197 204 ! = 2, if atoms i1 and i2 are separated by 3 covalent 198 205 ! bonds and their distance can change 199 206 ! = 1, for all other pairs 200 ! if abs(i2-i1) > mxconr, the atoms are assumed to be separated by 201 ! many bonds, and with no restriction on their distances. On a protein 207 ! if abs(i2-i1) > mxconr, the atoms are assumed to be separated by 208 ! many bonds, and with no restriction on their distances. On a protein 202 209 ! molecule made of natural amino acids, atoms with indices separated 203 210 ! by more than 35 can not be connected by three covalent bonds. … … 228 235 do m=1,nbdat(iat2) 229 236 iat3p=ibdat(m,iat2) 230 if (iat3p.ne.iat3) then 237 if (iat3p.ne.iat3) then 231 238 matcon(iat4-iat3p,iat3p)=2 ! 3 covalent bonds 232 239 matcon(iat3p-iat4,iat4)=2 ! … … 255 262 matcon(iat2-iat1,iat1)=0 256 263 matcon(iat1-iat2,iat2)=0 257 enddo 264 enddo 258 265 matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2 259 266 matcon(iCa(irsml1(iml))+1-iat1,iat1)=2 … … 266 273 ! Below: for certain residues, some atoms separated by 3 or more bonds 267 274 ! do not change distance. So, the connection matrix term for such pairs 268 ! should be zero. 275 ! should be zero. 269 276 270 277 do irs=irsml1(iml),irsml2(iml) 271 278 if (irs.eq.irsml1(iml)) then 272 279 iatoff=1 273 else 280 else 274 281 iatoff=0 275 282 endif … … 283 290 if ((mynm.eq.'pro')) then 284 291 prlvr=.true. ! residue i is a proline variant 285 else 292 else 286 293 prlvr=.false. 287 endif 288 if ((mynm.eq.'asn')) then 294 endif 295 if ((mynm.eq.'asn')) then 289 296 call set_local_constr(irs,5,9,2) 290 else if ((mynm.eq.'gln')) then 297 else if ((mynm.eq.'gln')) then 291 298 call set_local_constr(irs,8,12,5) 292 299 else if ((mynm.eq.'arg')) then … … 294 301 else if ((mynm.eq.'his')) then 295 302 call set_local_constr(irs,5,12,2) 296 else if (mynm.eq.'phe') then 303 else if (mynm.eq.'phe') then 297 304 call set_local_constr(irs,5,15,2) 298 305 else if (mynm.eq.'tyr') then … … 308 315 call set_local_constr(irs,5,19,2) 309 316 else if (prlvr) then 310 ! Proline. Many more distances are fixed because of the fixed 317 ! Proline. Many more distances are fixed because of the fixed 311 318 ! phi angle 312 319 call set_local_constr(irs,-3,11,-3) … … 347 354 do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml)) 348 355 do iat2=iat1+1,iatrs2(irsml2(iml)) 349 if (iat2-iat1.le.mxconr) then 350 if (matcon(iat2-iat1,iat1).eq.2) then 351 ilp=ilp+1 352 lcp1(ilp)=iat1 353 lcp2(ilp)=iat2 354 endif 356 if ((iat2-iat1.le.mxconr).and. 357 & matcon(iat2-iat1,iat1).eq.2) then 358 ilp=ilp+1 359 lcp1(ilp)=iat1 360 lcp2(ilp)=iat2 355 361 endif 356 362 enddo … … 358 364 359 365 ilpnd(iml)=ilp 360 if (iml.lt.ntlml) then 366 if (iml.lt.ntlml) then 361 367 ilpst(iml+1)=ilp+1 362 368 endif … … 369 375 enddo 370 376 enddo 377 371 378 print *,'finished initializing Lund force field' 372 379 end … … 374 381 integer function ihptype(iaa) 375 382 include 'INCL.H' 383 integer iaa, ityp 384 376 385 character mynm*4 377 386 mynm=seq(iaa) … … 383 392 & .or.(mynm.eq.'met').or.(mynm.eq.'pro')) then 384 393 ityp=2 385 else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp')) 394 else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp')) 386 395 & then 387 396 ityp=3 388 397 endif 389 398 ihptype=ityp 390 return 399 return 391 400 end 392 401 … … 394 403 include 'INCL.H' 395 404 include 'incl_lund.h' 405 double precision q1, q2, et, xij, yij, zij 406 407 integer i, iat1, irsd, j, iat2 408 396 409 dimension q1(2),q2(2) 397 410 data q1/-0.2,0.2/ 398 data q2/0.42,-0.42/ 411 data q2/0.42,-0.42/ 399 412 400 413 et=0.0 … … 412 425 return 413 426 end 414 ! Evaluates backbone backbone hydrogen bond strength for residues 427 ! Evaluates backbone backbone hydrogen bond strength for residues 415 428 ! i and j, taking the donor from residue i and acceptor from residue j 416 429 real*8 function ehbmmrs(i,j) 417 430 include 'INCL.H' 418 431 include 'incl_lund.h' 432 double precision rdon2, racc2, evlu 433 434 integer i, j 435 419 436 double precision r2,r4,r6,dx,dy,dz,ca,cb 420 437 integer d1,d2,a1,a2 … … 427 444 dz=zat(a1)-zat(d2) 428 445 r2=dx*dx+dy*dy+dz*dz 429 if (r2.gt.cthb2) then 446 if (r2.gt.cthb2) then 430 447 ! print *,'hbmm = 0 ',cthb2,r2,a1,a2,d1,d2 431 448 ! print *,'a1,a2,d1,d2,r2 = ',a1,a2,d1,d2,r2,sighb2,cthb … … 464 481 include 'INCL.H' 465 482 include 'incl_lund.h' 483 double precision ebiasrs, shbm1, shbm2, etmp, ehbmmrs, ehp, etmp1 484 double precision xij, yij, zij, r2, r6 485 486 integer istres, nml, indres, i, j, ihpi, ihptype, ihpj, i1, i2 487 integer iat1, iat2, iatt1, iatt2 488 466 489 character mynm*4 467 490 logical prlvr … … 485 508 & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then 486 509 prlvr=.true. ! residue i is a proline variant 487 else 510 else 488 511 prlvr=.false. 489 endif 512 endif 490 513 491 514 ! Bias, or local electrostatic term. Excluded from the list are … … 502 525 ! j run over the whole set of amino acids. 503 526 ! No terms for residue i, if it is a proline variant. 504 if (.not.prlvr) then 527 if (.not.prlvr) then 505 528 do j=istres,indres 506 if ((j.lt.(i-2)).or.(j.gt.(i+1))) then 529 if ((j.lt.(i-2)).or.(j.gt.(i+1))) then 507 530 shbm2=1.0 508 531 if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5 509 532 etmp=ehbmmrs(i,j) 510 533 eyhb=eyhb+shbm1*shbm2*etmp 511 endif 534 endif 512 535 enddo 513 536 endif 514 537 ! Hydrophobicity, only if residue i is hydrophobic to start with 515 538 ihpi=ihptype(i) 516 if (ihpi.ge.0) then 539 if (ihpi.ge.0) then 517 540 ! Unlike hydrogen bonds, the hydrophobicity potential is symmetric 518 541 ! in i and j. So, the loop for j runs from i+1 to the end. … … 526 549 ! print *, 'hydrophobicity contribution ',etmp,i,j 527 550 eysl=eysl+etmp 528 endif 551 endif 529 552 enddo 530 553 endif … … 535 558 ! Local pair or third-neighbour excluded volume 536 559 ! Numerically this is normally the term with largest positive 537 ! contribution to the energy in an equilibrated stystem. 560 ! contribution to the energy in an equilibrated stystem. 538 561 539 562 i1=ilpst(nml) … … 551 574 r2=(xij*xij + yij*yij + zij*zij) 552 575 if (r2.le.exvcut2) then 553 r6=sig2lcp(iatt1,iatt2)/r2 576 r6=sig2lcp(iatt1,iatt2)/r2 554 577 r6=r6*r6*r6 555 578 etmp1=(r6*r6+asalcp(iatt1,iatt2)+bsalcp(iatt1,iatt2)*r2) … … 560 583 etmp=etmp+etmp1 561 584 endif 562 ! print *,'pair : ',iat1,iat2,' contribution ',etmp1 585 ! print *,'pair : ',iat1,iat2,' contribution ',etmp1 563 586 ! print *,exvcut2,r2 564 587 enddo … … 579 602 include 'INCL.H' 580 603 include 'incl_lund.h' 604 double precision b2, a2, r2min, xij, yij, zij, dtmp, ssum 605 606 integer ihp1, ihp2, ni, i1, nj, i2, i, k, j, l 607 581 608 dimension r2min(12) 582 609 … … 608 635 ssum=0 609 636 do i=1,ni+nj 610 if (r2min(i).le.b2) then 611 if (r2min(i).lt.a2) then 637 if (r2min(i).le.b2) then 638 if (r2min(i).lt.a2) then 612 639 ssum=ssum+1 613 else 640 else 614 641 ssum=ssum+(b2-r2min(i))/(b2-a2) 615 642 endif … … 626 653 include 'incl_lund.h' 627 654 ! For multi-chain systems it makes little sense to split the calculation 628 ! of this term into an 'interaction part' and a contribution from 655 ! of this term into an 'interaction part' and a contribution from 629 656 ! individual molecules. So, normally this should always be called with 630 657 ! argument nml=0. Only for diagnostic reasons, you might want to find 631 ! the contribution from one molecule in a multi-chain system assuming 658 ! the contribution from one molecule in a multi-chain system assuming 632 659 ! there was no other molecule. 660 double precision xmin, ymin, zmin, xmax, ymax, zmax, sizex, sizey 661 double precision sizez, shiftx, shifty, shiftz, etmp, xij, yij 662 double precision zij, r2, r6, etmp1 663 664 integer nml, istat, indat, i, incell, ndx, ndy, ndz, nxy, ncell 665 integer nocccl, locccl, j, mx, my, mz, icellj, icell, jj, isort 666 integer icl, lcell, ix, iz, iy, nex, ney, nez, nsame, nngbr, jx 667 integer jy, jz, jcl, ii, ngbr, i1, iat1, i2, iat2, iatt1, iatt2 668 633 669 dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell) 634 670 dimension icell(mxat) 635 671 integer :: jx1, jy1, jz1 636 logical doit 637 638 if (nml.eq.0) then 672 673 if (nml.eq.0) then 639 674 istat=iatrs1(irsml1(1)) 640 675 indat=iatrs2(irsml2(ntlml)) 641 else 676 else 642 677 istat=iatrs1(irsml1(nml)) 643 678 indat=iatrs2(irsml2(nml)) … … 645 680 646 681 eyvw=0.0 647 ! The beginning part of this implementation is very similar to the 648 ! assignment of cells to the atoms during calculation of solvent 649 ! accessible surface area. So, much of that part is similar. But 682 ! The beginning part of this implementation is very similar to the 683 ! assignment of cells to the atoms during calculation of solvent 684 ! accessible surface area. So, much of that part is similar. But 650 685 ! unlike the accessible surface calculations, this term is symmetric 651 686 ! in any two participating atoms. So, the part after the assignment … … 724 759 icellj=mx+my*ndx+mz*nxy+1 725 760 icell(j)=icellj 726 if (icellj.gt.mxcell) then 761 if (icellj.gt.mxcell) then 727 762 print *,'exvlun> bad cell index ',icellj,' for atom ',j 728 763 stop 729 else 730 if (incell(icellj).eq.0) then 764 else 765 if (incell(icellj).eq.0) then 731 766 ! previously unoccupied cell 732 767 nocccl=nocccl+1 … … 756 791 ix=mod(lcell-1,ndx) 757 792 iz = (lcell - 1) / nxy 758 iy = ((lcell - 1) - iz * nxy) / ndx 793 iy = ((lcell - 1) - iz * nxy) / ndx 759 794 760 795 c print *,'icl=',icl,'absolute index of cell = ',lcell … … 767 802 nngbr=0 768 803 do jx=ix,nex 769 if (jx.ge.ndx) then 804 if (jx.ge.ndx) then 770 805 jx1=0 771 else 806 else 772 807 jx1=jx 773 808 endif 774 809 do jy=iy,ney 775 if (jy.ge.ndy) then 810 if (jy.ge.ndy) then 776 811 jy1=0 777 else 812 else 778 813 jy1=jy 779 endif 814 endif 780 815 do jz=iz,nez 781 if (jz.ge.ndz) then 816 if (jz.ge.ndz) then 782 817 jz1=0 783 else 818 else 784 819 jz1=jz 785 820 endif 786 821 jcl=jx1 + ndx*jy1 + nxy*jz1 + 1 787 ! write (*,*)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1822 ! write (logString, *)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1 788 823 do ii=incell(jcl)+1,incell(jcl+1) 789 824 ! count the total number of neighbours … … 799 834 enddo 800 835 ! A few more cells need to be searched, so that we cover 13 of the 26 801 ! neighbouring cells. 836 ! neighbouring cells. 802 837 ! 1 803 838 jx=ix+1 804 if (jx.ge.ndx) jx=0 839 if (jx.ge.ndx) jx=0 805 840 jy=iy 806 841 jz=iz-1 … … 883 918 r2=(xij*xij+yij*yij+zij*zij) 884 919 885 if (r2.le.exvcutg2) then 886 doit=.false. 887 if (abs(iat2-iat1).gt.mxconr ) then 888 doit=.true. 889 else if (matcon(iat2-iat1,iat1).eq.1) then 890 doit=.true. 891 endif 892 if (doit) then 920 if (r2.le.exvcutg2) then 921 if (abs(iat2-iat1).gt.mxconr.or. 922 & matcon(iat2-iat1,iat1).eq.1) then 893 923 iatt1=ityat(iat1) 894 924 iatt2=ityat(iat2) … … 931 961 include 'incl_lund.h' 932 962 963 double precision etmp, etmp1, xij, yij, zij, r2, r6 964 965 integer iat1, iat2, iatt1, iatt2 966 933 967 etmp=0.0 934 968 etmp1=0.0 … … 941 975 r2=(xij*xij+yij*yij+zij*zij) 942 976 943 if (r2.le.exvcutg2) then 977 if (r2.le.exvcutg2) then 944 978 if (abs(iat2-iat1).gt.mxconr.or. 945 979 & matcon(iat2-iat1,iat1).eq.1) then … … 951 985 & +bsaexv(iatt1,iatt2)*r2 952 986 etmp=etmp+etmp1 953 else 987 else 954 988 ! print *,'atoms ', iat1,' and ',iat2,' were close', 955 989 ! # 'but matcon is ',matcon(iat2-iat1,iat1) -
enyreg.f
r6650a56 r32289cd 21 21 22 22 23 double precision eny 24 25 integer i, nml, j 26 23 27 eny = 0.d0 24 28 -
enyshe.f
r6650a56 r32289cd 1 1 ! ************************************************************** 2 2 ! 3 ! This file contains the subroutines: enyshe 3 ! This file contains the subroutines: enyshe 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 20 20 ! 21 21 ! The function loops over all moving sets within the molecule. Within 22 ! this loop it loops over the van-der-Waals domains of each atom in the 22 ! this loop it loops over the van-der-Waals domains of each atom in the 23 23 ! moving set and finally over the atoms that belong to the 1-4 interaction 24 24 ! set. … … 28 28 29 29 ! If nml == 0 calculate the interaction between all pairs. 30 double precision e0, vr, cqi, xi, yi, zi, xij, yij, zij, rij2 31 double precision rij4, rij6, rij, sr, ep 32 33 integer nml, ntlvr, ifivr, i1s, io, iv, ia, it, ic, i2s, ims, i1 34 integer i2, i, ity, ivw, j, jty, i14 35 30 36 if (nml.eq.0) then 31 37 ntlvr = nvr … … 33 39 ntlvr=nvrml(nml) 34 40 endif 35 41 36 42 if (ntlvr.eq.0) then 37 43 write (*,'(a,i4)') … … 48 54 ifivr = ivrml1(1) 49 55 i1s = imsml1(ntlml) + nmsml(ntlml) 50 else 56 else 51 57 ! Index of first variable in molecule. 52 58 ifivr=ivrml1(nml) … … 54 60 i1s=imsml1(nml)+nmsml(nml) 55 61 endif 56 ! Loop over moving sets/variables in reverse order 57 do io=ifivr+ntlvr-1,ifivr,-1 62 ! Loop over moving sets/variables in reverse order 63 do io=ifivr+ntlvr-1,ifivr,-1 58 64 ! The array iorvr contains the variables in an "apropriate" order. 59 iv=iorvr(io) 65 iv=iorvr(io) 60 66 ! Index of the primary moving atom for the variable with index iv 61 ia=iatvr(iv) 67 ia=iatvr(iv) 62 68 ! Get the type of variable iv (valence length, valence angle, dihedral angle) 63 it=ityvr(iv) 69 it=ityvr(iv) 64 70 ! Class of variable iv's potential (Q: What are they) 65 ic=iclvr(iv) 71 ic=iclvr(iv) 66 72 ! If iv is a dihedral angle ... 67 if (it.eq.3) then 73 if (it.eq.3) then 68 74 ! Barrier height * 1/2 of the potential of iv. 69 75 e0=e0to(ic) 70 76 ! Calculate the periodic potential term. sgto is the sign of the barrier, rnto is 71 77 ! the periodicity and toat is torsion angle(?) associate with atom ia. 72 if (e0.ne.0.) 78 if (e0.ne.0.) 73 79 & eyvr=eyvr+e0*(1.0+sgto(ic)*cos(toat(ia)*rnto(ic))) 74 80 ! else if iv is a valence angle ... 75 elseif (it.eq.2) then 81 elseif (it.eq.2) then 76 82 ! vr is the valence angle of ia 77 83 vr=baat(ia) 78 84 ! else if iv is a valence length... 79 elseif (it.eq.1) then 85 elseif (it.eq.1) then 80 86 ! vr is the length of the valence bond 81 87 vr=blat(ia) … … 86 92 i2s=i1s-1 87 93 ! index of first moving set associated with iv 88 i1s=imsvr1(iv) 94 i1s=imsvr1(iv) 89 95 ! Loop over all moving sets starting from the one associated with vr to the end. 90 do ims=i1s,i2s 96 do ims=i1s,i2s 91 97 ! First atom of the current moving set 92 98 i1=latms1(ims) … … 94 100 i2=latms2(ims) 95 101 ! Loop over all atoms of the current moving set. 96 do i=i1,i2 102 do i=i1,i2 97 103 ! Atom class of current atom 98 104 ity=ityat(i) … … 104 110 zi=zat(i) 105 111 ! Loop over the atoms of the van der Waals domain belonging to atom i 106 do ivw=ivwat1(i),ivwat2(i) 107 ! Loop over the atoms of the van der Waals domain of the atoms of the 112 do ivw=ivwat1(i),ivwat2(i) 113 ! Loop over the atoms of the van der Waals domain of the atoms of the 108 114 ! van der Waals domain of atom i 109 115 ! Q: Which atoms are in these domains? 110 do j=lvwat1(ivw),lvwat2(ivw) 116 do j=lvwat1(ivw),lvwat2(ivw) 111 117 ! Atom type of partner 112 118 jty=ityat(j) … … 139 145 endif 140 146 141 enddo 142 enddo 143 147 enddo 148 enddo 149 144 150 ! Loop over 1-4 interaction partners 145 151 ! The interactions between atoms that are three bonds apart in the protein are 146 152 ! dominated by quantum mechanical effects. They are treated separately. 147 do i14=i14at1(i),i14at2(i) 153 do i14=i14at1(i),i14at2(i) 148 154 j=l14at(i14) 149 155 -
enyshe_p.f
r6650a56 r32289cd 1 1 ! ************************************************************** 2 2 ! 3 ! This file contains the subroutines: enyshe 3 ! This file contains the subroutines: enyshe 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 14 14 15 15 ! ............................................................................ 16 ! 16 ! 17 17 ! PURPOSE: Calculate internal energy of molecule 'nml' with ECEPP parameters 18 ! 18 ! 19 19 ! CALLS: none 20 ! 20 ! 21 21 ! The function loops over all moving sets within the molecule. Within 22 ! this loop it loops over the van-der-Waals domains of each atom in the 22 ! this loop it loops over the van-der-Waals domains of each atom in the 23 23 ! moving set and finally over the atoms that belong to the 1-4 interaction 24 24 ! set. … … 31 31 32 32 ! If nml == 0 calculate the interaction between all pairs. 33 double precision eysmsum, teysm, teyel, teyvw, teyhb, teyvr 34 double precision startwtime, e0, vr, cqi, xi, yi, zi, xij, yij 35 double precision zij, rij2, rij4, rij6, rij, sr, ep, endwtime 36 37 integer nml, ntlvr, ifivr, i1s, iend, istart, loopcounter, io, iv 38 integer ia, it, ic, i2s, ims, i1, i2, i, ity, ivw, j, jty, i14 39 integer ierror 40 33 41 if (nml.eq.0) then 34 42 ntlvr = nvr … … 36 44 ntlvr=nvrml(nml) 37 45 endif 38 46 39 47 if (ntlvr.eq.0) then 40 48 write (*,'(a,i4)') … … 56 64 ifivr = ivrml1(1) 57 65 i1s = imsml1(ntlml) + nmsml(ntlml) 58 else 66 else 59 67 ! Index of first variable in molecule. 60 68 ifivr=ivrml1(nml) … … 62 70 i1s=imsml1(nml)+nmsml(nml) 63 71 endif 64 ! Loop over variables in reverse order 72 ! Loop over variables in reverse order 65 73 ! This is the first loop to parallize. We'll just split the moving sets 66 74 ! over the number of available processors and sum the energy up in the end. … … 72 80 startwtime = MPI_Wtime() 73 81 loopcounter = 0 74 ! do io=ifivr+ntlvr-1,ifivr,-1 75 do io = workPerProcessor(nml, myrank) - 1, 82 ! do io=ifivr+ntlvr-1,ifivr,-1 83 do io = workPerProcessor(nml, myrank) - 1, 76 84 & workPerProcessor(nml, myrank+1), -1 77 85 if (io.lt.istart) then … … 79 87 endif 80 88 ! The array iorvr contains the variables in an "apropriate" order. 81 iv=iorvr(io) 89 iv=iorvr(io) 82 90 ! Index of the primary moving atom for the variable with index iv 83 ia=iatvr(iv) 91 ia=iatvr(iv) 84 92 ! Get the type of variable iv (valence length, valence angle, dihedral angle) 85 it=ityvr(iv) 93 it=ityvr(iv) 86 94 ! Class of variable iv's potential (Q: What are they) 87 ic=iclvr(iv) 95 ic=iclvr(iv) 88 96 ! If iv is a dihedral angle ... 89 if (it.eq.3) then 97 if (it.eq.3) then 90 98 ! Barrier height * 1/2 of the potential of iv. 91 99 e0=e0to(ic) 92 100 ! Calculate the periodic potential term. sgto is the sign of the barrier, rnto is 93 101 ! the periodicity and toat is torsion angle(?) associate with atom ia. 94 if (e0.ne.0.) 102 if (e0.ne.0.) 95 103 & teyvr=teyvr+e0*(1.0+sgto(ic)*cos(toat(ia)*rnto(ic))) 96 104 ! else if iv is a valence angle ... 97 elseif (it.eq.2) then 105 elseif (it.eq.2) then 98 106 ! vr is the valence angle of ia 99 107 vr=baat(ia) 100 108 ! else if iv is a valence length... 101 elseif (it.eq.1) then 109 elseif (it.eq.1) then 102 110 ! vr is the length of the valence bond 103 111 vr=blat(ia) … … 108 116 i2s=i1s-1 109 117 ! index of first moving set associated with iv 110 i1s=imsvr1(iv) 118 i1s=imsvr1(iv) 111 119 ! Loop over all moving sets starting from the one associated with vr to the end. 112 do ims=i1s,i2s 120 do ims=i1s,i2s 113 121 ! First atom of the current moving set 114 122 i1=latms1(ims) … … 116 124 i2=latms2(ims) 117 125 ! Loop over all atoms of the current moving set. 118 do i=i1,i2 126 do i=i1,i2 119 127 ! Atom class of current atom 120 128 ity=ityat(i) … … 126 134 zi=zat(i) 127 135 ! Loop over the atoms of the van der Waals domain belonging to atom i 128 do ivw=ivwat1(i),ivwat2(i) 129 ! Loop over the atoms of the van der Waals domain of the atoms of the 136 do ivw=ivwat1(i),ivwat2(i) 137 ! Loop over the atoms of the van der Waals domain of the atoms of the 130 138 ! van der Waals domain of atom i 131 139 ! Q: Which atoms are in these domains? 132 do j=lvwat1(ivw),lvwat2(ivw) 140 do j=lvwat1(ivw),lvwat2(ivw) 133 141 134 142 loopcounter = loopcounter + 1 … … 163 171 endif 164 172 165 enddo 166 enddo 167 173 enddo 174 enddo 175 168 176 ! Loop over 1-4 interaction partners 169 177 ! The interactions between atoms that are three bonds apart in the protein are 170 178 ! dominated by quantum mechanical effects. They are treated separately. 171 do i14=i14at1(i),i14at2(i) 179 do i14=i14at1(i),i14at2(i) 172 180 loopcounter = loopcounter + 1 173 181 j=l14at(i14) … … 210 218 211 219 ! Collect energies from all nodes and sum them up 212 call MPI_ALLREDUCE(teysm, eysmsum, 1, MPI_DOUBLE_PRECISION, 220 call MPI_ALLREDUCE(teysm, eysmsum, 1, MPI_DOUBLE_PRECISION, 213 221 & MPI_SUM, my_mpi_comm, ierror) 214 222 call MPI_ALLREDUCE(teyel, eyel, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 215 223 & my_mpi_comm, ierror) 216 call MPI_ALLREDUCE(teyvw, eyvw, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 217 & my_mpi_comm, ierror) 218 call MPI_ALLREDUCE(teyhb, eyhb, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 219 & my_mpi_comm, ierror) 220 call MPI_ALLREDUCE(teyvr, eyvr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 224 call MPI_ALLREDUCE(teyvw, eyvw, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 225 & my_mpi_comm, ierror) 226 call MPI_ALLREDUCE(teyhb, eyhb, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 227 & my_mpi_comm, ierror) 228 call MPI_ALLREDUCE(teyvr, eyvr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 221 229 & my_mpi_comm, ierror) 222 230 -
enysol.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 10 10 ! ************************************************************** 11 11 12 12 13 13 real*8 function enysol(nmol) 14 14 … … 17 17 ! 18 18 ! Double Cubic Lattice algorithm for calculating the 19 ! solvation energy of proteins using 19 ! solvation energy of proteins using 20 20 ! solvent accessible area method. 21 21 ! … … 24 24 ! 25 25 ! ------------------------------------------------------------- 26 ! TODO: Check the solvent energy for multiple molecules 26 ! TODO: Check the solvent energy for multiple molecules 27 27 ! arguments 28 28 integer nmol 29 29 30 30 ! functions 31 31 integer nursat … … 41 41 double precision sizes, trad, zmin, xmax, xmin, ymax, zmax 42 42 double precision sdr, sdd, volume 43 43 44 44 dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat) 45 45 dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox) … … 47 47 48 48 ! common/ressurf/surfres(mxrs) 49 49 50 50 eyslh = 0.0 51 51 eyslp = 0.0 … … 61 61 do i=nrslow,nrshi 62 62 surfres(i) = 0.0d0 63 end do 63 end do 64 64 65 65 numat= nup - nlow + 1 … … 68 68 inbox(i)=0 69 69 end do 70 70 71 71 asa=0.0d0 72 72 vdvol=0.0d0 … … 135 135 stop 136 136 end if 137 137 138 138 ! Let us shift the borders to home all boxes 139 139 … … 175 175 inbox(i+1)=inbox(i+1)+inbox(i) 176 176 end do 177 178 177 178 179 179 ! Sorting the atoms by the their box numbers 180 180 … … 184 184 indsort(jj)=i 185 185 inbox(j)=jj-1 186 end do 187 186 end do 187 188 188 ! Getting started 189 189 … … 205 205 ney=min(iy+1,ndy-1) 206 206 nez=min(iz+1,ndz-1) 207 207 208 208 ! Atoms in the boxes around 209 209 … … 221 221 end do 222 222 end do 223 end do 224 223 end do 224 225 225 do ia=inbox(ibox)+1,inbox(ibox+1) 226 226 jbi=indsort(ia) … … 274 274 end do 275 275 99 continue 276 276 277 277 if(ik.gt.nnei)then 278 278 surfc(il)=.true. … … 316 316 surfres(jres) = surfres(jres) + area 317 317 end if 318 end do 318 end do 319 319 end if 320 320 end do … … 341 341 character lin*80 342 342 integer i 343 ! Skipping comment lines, which begin with '!' 343 ! Skipping comment lines, which begin with '!' 344 344 read(20,'(a)') lin 345 345 do while(lin(1:1).eq.'!') … … 352 352 ! write(*,'(a,i5)') 'the number of points---->',npnt 353 353 354 ! Read the surface points 354 ! Read the surface points 355 355 356 356 do i=1,npnt … … 358 358 ! write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3) 359 359 end do 360 360 361 361 return 362 362 363 363 end 364 364 -
enysol_p.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 10 10 ! ************************************************************** 11 11 12 12 13 13 real*8 function enysol(nmol) 14 14 … … 19 19 ! 20 20 ! Double Cubic Lattice algorithm for calculating the 21 ! solvation energy of proteins using 21 ! solvation energy of proteins using 22 22 ! solvent accessible area method. 23 23 ! … … 26 26 ! 27 27 ! ------------------------------------------------------------- 28 ! TODO: Check the solvent energy for multiple molecules 28 ! TODO: Check the solvent energy for multiple molecules 29 double precision startwtime, avr_x, avr_y, avr_z, xmin, ymin, zmin 30 double precision xmax, ymax, zmax, rmax, diamax, sizex, sizey 31 double precision sizez, shiftx, shifty, shiftz, boxpp, trad, dx 32 double precision dy, dz, dd, akrad, dr, xyz, radb, radb2, sdd, sdr 33 double precision area, volume, eyslsum, endwtime 34 35 integer nmol, nrslow, nrshi, nlow, nup, i, numat, inbox, j, ndx 36 integer ndy, ndz, nqxy, ncbox, mx, my, mz, nboxj, numbox, jj 37 integer indsort, iboxmin, iboxmax, ibox, iz, iy, ix, lbn, nsx, nsy 38 integer nsz, nex, ney, nez, jcnt, jz, jy, jx, jbox, ii, look, ia 39 integer jbi, nnei, ib, jtk, il, lst, ilk, ik, icount, jres, nursat 40 integer ierror, nhx, mhx, nbt, mbt 41 29 42 dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat) 30 43 dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox) … … 33 46 34 47 ! common/ressurf/surfres(mxrs) 35 real*8 tsurfres(mxrs) 48 real*8 tsurfres(mxrs) 36 49 startwtime = MPI_Wtime() 37 50 root = 0 … … 50 63 do i=nrslow,nrshi 51 64 surfres(i) = 0.0d0 52 end do 65 end do 53 66 54 67 numat= nup - nlow + 1 … … 57 70 inbox(i)=0 58 71 end do 59 72 60 73 asa=0.0d0 61 74 vdvol=0.0d0 … … 124 137 stop 125 138 end if 126 139 127 140 ! Let us shift the borders to home all boxes 128 141 … … 164 177 inbox(i+1)=inbox(i+1)+inbox(i) 165 178 end do 166 167 179 180 168 181 ! Sorting the atoms by the their box numbers 169 182 … … 173 186 indsort(jj)=i 174 187 inbox(j)=jj-1 175 end do 176 188 end do 189 177 190 ! Getting started 178 191 ! We have to loop over ncbox boxes and have no processors available … … 202 215 ney=min(iy+1,ndy-1) 203 216 nez=min(iz+1,ndz-1) 204 217 205 218 ! Atoms in the boxes around 206 219 … … 218 231 end do 219 232 end do 220 end do 221 233 end do 234 222 235 do ia=inbox(ibox)+1,inbox(ibox+1) 223 236 jbi=indsort(ia) … … 271 284 end do 272 285 99 continue 273 286 274 287 if(ik.gt.nnei)then 275 288 surfc(il)=.true. … … 313 326 surfres(jres) = surfres(jres) + area 314 327 end if 315 end do 328 end do 316 329 end if 317 330 ! end do 318 331 ! end do 319 332 end do 320 call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION, 333 call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION, 321 334 & MPI_SUM,my_mpi_comm, ierror) 322 335 ! write(*,*) 'enysol>', myrank, eysl, eyslsum 323 336 tsurfres = surfres 324 call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION, 337 call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION, 325 338 & MPI_SUM,my_mpi_comm, ierror) 326 339 eysl = eyslsum 327 340 328 341 endwtime = MPI_Wtime() 329 342 if (myrank.le.-1) then … … 349 362 subroutine tessel 350 363 include 'INCL.H' 364 integer i 365 351 366 character lin*80 352 367 353 ! Skipping comment lines, which begin with '!' 368 ! Skipping comment lines, which begin with '!' 354 369 355 370 read(20,'(a)') lin … … 363 378 ! write(*,'(a,i5)') 'the number of points---->',npnt 364 379 365 ! Read the surface points 380 ! Read the surface points 366 381 367 382 do i=1,npnt 368 383 read(20,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3) 369 384 370 385 ! write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3) 371 386 end do 372 387 373 388 return 374 389 375 390 end 376 391 -
esolan.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 14 14 ! ----------------------------------------------------------------- 15 15 ! 16 ! Calculates the solvation energy of the protein with 17 ! solavtion parameters model E=\sum\sigma_iA_i. 16 ! Calculates the solvation energy of the protein with 17 ! solavtion parameters model E=\sum\sigma_iA_i. 18 18 ! The solvent accessible surface area per atom 19 ! and its gradients with respect to the Cartesian coordinates of 19 ! and its gradients with respect to the Cartesian coordinates of 20 20 ! the atoms are calculated using the projection method described in 21 21 ! … … 25 25 ! 26 26 ! INPUT: nmol - the order number of the protein chain. 27 ! The atomic coordinates, radii and solvation parameters 27 ! The atomic coordinates, radii and solvation parameters 28 28 ! are taken from the global arrays of SMMP 29 29 ! 30 ! OUTPUT: global array gradan(mxat,3) contains the gradients of 30 ! OUTPUT: global array gradan(mxat,3) contains the gradients of 31 31 ! solvation energy per atom 32 32 ! local array as(mxat) contains accessible surface area per atom … … 42 42 43 43 include 'INCL.H' 44 45 44 45 46 46 integer nmol, ks0, ks2 47 47 Cf2py intent(in) nmol 48 48 49 49 parameter (ks0=mxat*(mxat+1)) 50 50 parameter (ks2=mxat+mxat) … … 53 53 54 54 integer neib, neibp, ivrx, iv, ial, i, ij, j, iii, idi,ivr, ifu, 55 & ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll, 55 & ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll, 56 56 & nb, numat, nyx 57 57 real*8 vertex, ax, pol, as, ayx, ayx1, probe, dd, ddat, dadx, gp, 58 58 & dalp,dbet, daalp, dabet, vrx, dv, dx, dy, dz, dt, di, 59 & dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold, 59 & dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold, 60 60 & energy, sss, a1, a, ab1, a2, aa, ac2, ac1, ab2, am2, 61 61 & am, aia, ak, alp, am1, ap1, amibe, amabe, amaal, amial, … … 63 63 & ba, bcd, bc, bet, c1, c, c2, cc, sfcg, caba, cg, cfsg, 64 64 & cfcg, daba, ct, cs, d1, d2, dbe, dal, dcotbma, dcbet, 65 & dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp, 66 & dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22, 65 & dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp, 66 & dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22, 67 67 & r22p, r2, pom, prob, r, r42p, ratom, rr, s, sfsg, sf, 68 & sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv, 68 & sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv, 69 69 & zz, zv 70 70 dimension neib(0:mxat),neibp(0:mxat),ivrx(ks0) … … 78 78 integer ta2(0:mxat),ta3(0:mxat),fullarc(0:mxat),al(0:ks2) 79 79 real*8 neibor(mxat,4) 80 80 81 81 real*8, dimension(:, :, :), allocatable :: grad 82 82 83 83 allocate(grad(mxat, mxat, 3)) 84 85 84 85 86 86 ! open(3,file='solvation.dat')! Output file with solvation data 87 87 … … 90 90 if (nmol.eq.0) then 91 91 numat=iatrs2(irsml2(ntlml))-iatrs1(irsml1(1))+1 ! Number of atoms 92 else 92 else 93 93 numat=iatrs2(irsml2(nmol))-iatrs1(irsml1(nmol))+1 ! Number of atoms 94 94 endif … … 96 96 do i=1,numat 97 97 do j=1,3 98 gradan(i,j)=0.0d0 98 gradan(i,j)=0.0d0 99 99 end do 100 100 end do … … 102 102 do i=1,numat 103 103 xold(i)=xat(i)! Redefine the coordinates just for safety 104 yold(i)=yat(i)! 105 zold(i)=zat(i)! 106 pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main' 104 yold(i)=yat(i)! 105 zold(i)=zat(i)! 106 pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main' 107 107 As(i)=pi4*pol(i)! Initially the whole surface of the atom is accessible. 108 108 end do 109 109 110 110 do iii=1,numat 111 111 do jjj=1,numat … … 120 120 !*************************** 121 121 122 do 1 i=1,numat ! Lop over atoms "i" 122 do 1 i=1,numat ! Lop over atoms "i" 123 123 ! ---------------------------------------------------- 124 125 R=rvdw(i) 124 125 R=rvdw(i) 126 126 jj=0 127 127 do j=1,numat !Find the neighbours of "i" … … 132 132 write(*,*)'ERROR in data: centres of two atoms coincide!' 133 133 write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i), 134 & xold(j),yold(j),zold(j),rvdw(j) 134 & xold(j),yold(j),zold(j),rvdw(j) 135 135 stop 'Centres of atoms coincide!' 136 136 endif … … 146 146 end if 147 147 end do!Neighbours 148 148 149 149 neib(0)=jj ! The number of neighbors of "i" 150 150 … … 197 197 ! 198 198 uga=grnd() ! Random \gamma angle 199 ufi=grnd() ! Random \pi angle 199 ufi=grnd() ! Random \pi angle 200 200 uga=pi*uga/2 201 201 ufi=pi*ufi/2 … … 219 219 end do 220 220 end do 221 ! 222 221 ! 222 223 223 if(jjj.ne.0) then ! Rotation is necessary 224 224 sfsg=sf*sg … … 243 243 do jj=1,neib(0) 244 244 j=neib(jj) 245 neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+ 245 neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+ 246 246 & (ddat(j,4)-R)**2-pol(j) ! a 247 247 neibor(j,2)=-pom*ddat(j,2) ! b … … 259 259 do while(k.gt.1) ! B 260 260 k=k-1 261 ! Analyse mutual disposition of every pair of neighbours 261 ! Analyse mutual disposition of every pair of neighbours 262 262 do 13 L=neib(0),k+1,-1 ! A 263 263 … … 283 283 elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le. ! a01 02 284 284 & -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then 285 ! The circle neib(k) encloses circle neib(L) and the later is discarded 285 ! The circle neib(k) encloses circle neib(L) and the later is discarded 286 286 neib(0)=neib(0)-1 287 287 do j=L,neib(0) … … 326 326 elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge. ! a03 02 327 327 & dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then 328 ! Don't exclude neib(k) 328 ! Don't exclude neib(k) 329 329 neib(0)=neib(0)-1 330 330 do j=k,neib(0) … … 370 370 elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge. ! a07 02 371 371 & dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then 372 ! discard the circle neib(L) 372 ! discard the circle neib(L) 373 373 neib(0)=neib(0)-1 374 374 do j=L,neib(0) … … 385 385 p2=(b1-b2)*pom 386 386 !-- Assign t and s coordinates and the order number of circles 387 vertex(iv,1)=(t+p1)/am 387 vertex(iv,1)=(t+p1)/am 388 388 vertex(iv,2)=(s-p2)/am 389 389 vertex(iv,3)=neib(k) … … 394 394 vertex(iv,3)=neib(k) 395 395 vertex(iv,4)=neib(L) 396 !-- 396 !-- 397 397 endif ! 10 398 398 … … 491 491 dbe=amabe-amibe 492 492 if(dal.lt.pi) then 493 As(i)=R22*(pi-dal) 493 As(i)=R22*(pi-dal) 494 494 elseif(dbe.lt.pi) then 495 495 As(i)=R22*(pi-dbe) … … 510 510 enddo 511 511 do j=1,neib(0) 512 ! "iv" is the number of intersection point. 512 ! "iv" is the number of intersection point. 513 513 do k=1,iv 514 514 if(neibp(j).eq.vertex(k,3).or.neibp(j).eq.vertex(k,4)) goto 30 … … 532 532 do k=1,ifu 533 533 a=neibor(fullarc(k),1) 534 if(a.lt.0.d0) then 534 if(a.lt.0.d0) then 535 535 aia=-1.d0 536 536 else … … 579 579 ak=a*a 580 580 daba=dabs(a) 581 if(a.lt.0d0) then 581 if(a.lt.0d0) then 582 582 aia=-1d0 583 583 else … … 679 679 do j=1,nyx 680 680 681 ! 681 ! 682 682 ! Escape 'bad' intersections. 683 683 if(dabs(ayx1(j+1)-ayx1(j)).lt.1d-8) goto 40 … … 686 686 ap1=ct+rr*dcos(prob) ! The middle point of the arc. 687 687 ap2=cs+rr*dsin(prob) 688 ! Verify if the middle point belongs to a covered part. 688 ! Verify if the middle point belongs to a covered part. 689 689 ! If yes then omit. 690 690 do ij=1,ial … … 754 754 dabet(3)=(a*ac2-ak*c2+wv*daba*a2*dsbet)/am2 755 755 dabet(4)=(2d0*ak*a2)/am2 756 dv(2)=R42*b*vv1 757 dv(3)=R42*c*vv1 758 dv(4)=(d-R42*a)*vv1 759 dv(1)=-dv(4)*R42 760 dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2 761 dx(2)=-(d+R42*a)*dv(2)*vv2 762 dx(3)=-(d+R42*a)*dv(3)*vv2 763 dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2 756 dv(2)=R42*b*vv1 757 dv(3)=R42*c*vv1 758 dv(4)=(d-R42*a)*vv1 759 dv(1)=-dv(4)*R42 760 dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2 761 dx(2)=-(d+R42*a)*dv(2)*vv2 762 dx(3)=-(d+R42*a)*dv(3)*vv2 763 dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2 764 764 dy(1)=(-2d0*d+4d0*R42*a)/daba*vv1-(b*b+c*c-2d0*a*d+2d0* 765 765 & R42*ak)*(aia*vv+daba*dv(1))/ak*vv2 … … 769 769 & R42*ak)*aia*dv(3)/a*vv2 770 770 dy(4)=-2.d0*a/daba*vv1-(b*b+c*c-2.d0*a*d+2.d0* 771 & R42*ak)*aia*dv(4)/a*vv2 771 & R42*ak)*aia*dv(4)/a*vv2 772 772 dz(1)=-2d0*d/wv/a*vv1-wv*(vv+a*dv(1))/ak*vv2 773 773 dz(2)=b/wv/a*vv1-wv*dv(2)/a*vv2 … … 820 820 dii(4,2)=dii(1,2)*R42 821 821 dii(4,3)=2d0*(ddat(jk,4)+R)*R42 822 822 823 823 do idi=1,3 824 824 grad(i,jk,idi)=grad(i,jk,idi)+ … … 826 826 & di(3)*dii(3,idi)+di(4)*dii(4,idi) 827 827 828 enddo 828 enddo 829 829 do idi=1,4 830 830 dta(idi)=5d-1*daalp(idi)*(((yv-zv*((-b*dsbetpalp … … 838 838 di1(idi)=R2*(aia*daalp(idi)-4d0*xv*dta(idi)/(4d0+tt*tt)) 839 839 di2(idi)=R2*(-aia*dabet(idi)-4d0*xv*dtb(idi)/(4d0+tt*tt)) 840 enddo 840 enddo 841 841 842 842 dii(1,1)=2d0*ddat(kk,2) … … 866 866 40 continue 867 867 enddo 868 868 869 869 else ! analyse the case with lines (not necessary after rotation). 870 870 … … 920 920 enddo 921 921 probe(nyx+1)=ayx1(nyx)+1d0 922 922 923 923 924 924 bc=b*b+c*c … … 960 960 22 continue 961 961 endif 962 enddo 962 enddo 963 963 11 continue 964 964 … … 1017 1017 111 continue 1018 1018 1019 ! write(3,*)' No Area sigma Enrg gradx grady', 1019 ! write(3,*)' No Area sigma Enrg gradx grady', 1020 1020 ! & ' gradz Rad Atom' 1021 1021 ! write(3,*) 1022 1022 1023 1023 1024 1024 do i=1,numat 1025 1025 enr=As(i)*sigma(i) … … 1038 1038 1039 1039 esolan=sss 1040 1040 1041 1041 deallocate(grad) 1042 1042 ! close(3) … … 1046 1046 200 format(2i4,3f16.6) 1047 1047 700 format(i5,7f8.3,5x,a4) 1048 1048 1049 1049 end -
eyabgn.f
r6650a56 r32289cd 5 5 ! Jan H. Meinke, Sandipan Mohanty 6 6 ! 7 ! Corrections to ECEPP energy terms due to R. A. Abagyan et al. 8 ! 7 ! Corrections to ECEPP energy terms due to R. A. Abagyan et al. 8 ! 9 9 ! Two terms are calculated: eyrccr and eyentr, representing respectively 10 ! c a term to slightly shift the backbone dihedral angle preferences in 10 ! c a term to slightly shift the backbone dihedral angle preferences in 11 11 ! the ECEPP potential slightly away from the helix region, and another 12 ! term to estimate the side-chain entropy from a given configuration. 12 ! term to estimate the side-chain entropy from a given configuration. 13 13 ! 14 14 ! … … 16 16 real*8 function eyrccr(nml) 17 17 include 'INCL.H' 18 double precision et, rsscl 19 20 integer nml, istres, indres, i, ipsi, in, ica, ic, mlvr, iphi 21 18 22 dimension iN(mxrs),iCa(mxrs),iC(mxrs),mlvr(mxvr) 19 23 dimension iphi(mxrs),ipsi(mxrs) … … 34 38 call tolost(mynm) 35 39 if ((mynm.eq.'val').or.(mynm.eq.'ile').or. 36 & (mynm.eq.'thr')) then 40 & (mynm.eq.'thr')) then 37 41 rsscl=1.0 38 else 42 else 39 43 rsscl=0.5 40 44 endif … … 48 52 eyrccr=et 49 53 return 50 end 54 end 51 55 52 56 subroutine init_abgn … … 54 58 ! dimension rsstrg(mxrs) 55 59 ! common /abgncor/rsstrg 60 double precision xarea, estrg 61 62 integer i, istres, indres, imytyp, j 63 56 64 dimension xarea(nrsty),estrg(nrsty) 57 65 character mynm*4 … … 102 110 ! print *,'comparing ',mynm,' with ',rsnmcd(j),imytyp 103 111 enddo 104 if (imytyp.eq.0) then 112 if (imytyp.eq.0) then 105 113 print *, 'Unknown residue type ',seq(i) 106 114 print *, 'Abagyan term strength set to 0' 107 115 rsstrg(i)=0.0 108 else 116 else 109 117 rsstrg(i)=estrg(imytyp)/xarea(imytyp) 110 118 endif 111 119 ! print *,'residue ',i,seq(i),' type ',imytyp 112 ! print *, 'strength for residue ',i,seq(i),' is ',rsstrg(i) 120 ! print *, 'strength for residue ',i,seq(i),' is ',rsstrg(i) 113 121 enddo 114 122 print *, 'initialized Abagyan corrections to ECEPP force field' 115 123 end 116 124 117 125 real*8 function eyentr(nml) 118 126 include 'INCL.H' … … 121 129 ! common/ressurf/surfres(mxrs) 122 130 ! common/bet/beta 131 double precision eentr, aars, strh 132 133 integer nml, istres, indres, i 134 123 135 eentr=0 124 136 if (nml.eq.0) then … … 136 148 ! The maximal burial entropies were estimated at temperature 300k 137 149 ! The values in the array estrg are k_B * T (=300k) * Entropy 138 ! Presently we need it at temperature 1/beta, so we need to 150 ! Presently we need it at temperature 1/beta, so we need to 139 151 ! multiply the strengths in estrg with (1/beta)/(300 kelvin) 140 152 ! 300 kelvin is approximately 0.59576607 kcal/mol. … … 149 161 return 150 162 end 151 152 real*8 function eyabgn(nml) 163 164 real*8 function eyabgn(nml) 153 165 include 'INCL.H' 166 double precision eyrccr, eyentr 167 168 integer nml 169 154 170 eyabgn=eyrccr(nml)+eyentr(nml) 155 171 ! print *,'Abagyan term = ',eyabgn -
gradient.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 21 21 include 'INCL.H' 22 22 23 24 double precision esm 25 26 integer i, ivr1, ivr2, j 23 27 24 28 esm = 0.d0 -
init_energy.f
r6650a56 r32289cd 5 5 ! 6 6 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 7 ! Shura Hayryan, Chin-Ku 7 ! Shura Hayryan, Chin-Ku 8 8 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 9 9 ! Jan H. Meinke, Sandipan Mohanty … … 17 17 ! PURPOSE: initialize energy parameters 18 18 ! 0 => ECEPP2 or ECEPP3 depending on the value of sh2 19 ! 1 => FLEX 19 ! 1 => FLEX 20 20 ! 2 => Lund force field 21 21 ! 3 => ECEPP with Abagyan corrections … … 29 29 include 'INCL.H' 30 30 31 integer ll, iendst, its 32 31 33 character libdir*(*),tesfil*80 32 34 33 35 if (ientyp.eq.1) then 34 36 flex = .true. 35 else 37 else 36 38 flex = .false. 37 39 end if 38 40 39 41 lunlib=10 40 42 ll=iendst(libdir) … … 56 58 endif 57 59 58 if (ientyp .eq. 2) then 60 if (ientyp .eq. 2) then 59 61 reslib=libdir(1:ll)//'lib.lun' 60 62 endif … … 78 80 79 81 else 80 82 81 83 if (itysol.ne.0) then 82 84 write(*,'(a)') ' init_energy> undefined solvent type !' … … 110 112 111 113 include 'INCL.H' 114 115 double precision hbc, hba, ri, ai, aei, aic, a, aj, c, rij 116 117 integer i, j, iac, ido, jac, jdo 112 118 113 119 dimension hbc(mxhbdo,mxhbac),hba(mxhbdo,mxhbac) … … 308 314 309 315 include 'INCL.H' 316 317 integer i 310 318 311 319 ! Atom types ------------------------------------------------------------ … … 517 525 & 'ile ','leu ','lys ','lys+','met ','phe ','cpro','pro ','cpru', 518 526 & 'prou','pron','pro+','ser ','thr ','trp ','tyr ','val ' / 519 527 520 528 data (onltcd(i),i=1,nrsty)/ ! One-letter codes for amino acid types 521 529 & 'A', 'R', 'R', 'N', 'D', 'D', 'C', 'C', 'Q', … … 527 535 ! coefficients for their solvation free energy (kcal/molxA**2) 528 536 529 ! Method: 530 531 ! itysol=1 : OONS --> T.Ooi, et al, 537 ! Method: 538 539 ! itysol=1 : OONS --> T.Ooi, et al, 532 540 ! Proc. Natl. Acad. Sci. USA 8 (1987) 3086-3090. 533 ! itysol=2 : JRF --> J.Vila, et al, 541 ! itysol=2 : JRF --> J.Vila, et al, 534 542 ! PROTEINS: Struct Funct Genet 10(1991) 199-218. 535 ! itysol=3 : WE92 --> L.Wesson, D.Eisenberg, 543 ! itysol=3 : WE92 --> L.Wesson, D.Eisenberg, 536 544 ! Protein Science 1 (1992) 227-235. 537 545 ! itysol=4 : SCH1 --> D.Eisenberg, et al, … … 539 547 ! itysol=5 : SCH2 --> A.H.Juffer, et al, 540 548 ! Proteine Science 4 (1995) 2499-2509. 541 ! itysol=6 : SCH3 --> L.Wesson, D.Eisenberg, 549 ! itysol=6 : SCH3 --> L.Wesson, D.Eisenberg, 542 550 ! Protein Science 1 (1992) 227-235. 543 551 ! itysol=7 : SCH4 --> C.A. Schiffer, et al, 544 552 ! Mol. Simul. 10(1993) 121-149. 545 ! itysol=8 : EM86 --> D.Eisenberg, A.D. Mclachlan, 553 ! itysol=8 : EM86 --> D.Eisenberg, A.D. Mclachlan, 546 554 ! Nature 319 (1986) 199-203. 547 555 ! itysol=9 : BM --> B. Freyberg, et al, -
init_molecule.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 16 16 ! PURPOSE: construct starting structure of molecule(s) 17 17 ! 18 ! iabin = 1 : ab Initio using sequence & 18 ! iabin = 1 : ab Initio using sequence & 19 19 ! variables given in input files 20 20 ! iabin != 1 : sequence, variable information … … 32 32 include 'INCP.H' 33 33 34 integer iabin, iendst, ntl, i, j, l, it, ier, ir, nursvr, i1, i2 35 integer its 36 34 37 Cf2py character*80 optional, intent(in) :: seqfile = ' ' 35 Cf2py character*80 optional, intent(in) :: varfile = ' ' 36 38 Cf2py character*80 optional, intent(in) :: varfile = ' ' 39 37 40 character grpn*4,grpc*4 38 character navr*3, nars*4 41 character navr*3, nars*4 39 42 character seqfile*80, varfile*80 40 43 integer ontlml … … 44 47 readFromStdin = .false. 45 48 46 write ( *,*) 'init_molecule: Solvent: ', itysol47 if (iabin.eq.1) then 48 49 write (logString, *) 'init_molecule: Solvent: ', itysol 50 if (iabin.eq.1) then 51 49 52 ! ----------------------------------------- get sequence for molecule(s) 50 53 lunseq=11 … … 53 56 endif 54 57 if (iendst(seqfile).le.1.or.seqfile.eq.' ') then 55 1 write ( *,'(/,a,$)') ' file with SEQUENCE:'58 1 write (logString, '(/,a,$)') ' file with SEQUENCE:' 56 59 seqfil=' ' 57 60 read (*,'(a)',err=1) seqfil 58 61 readFromStdin = .true. 59 else 62 else 60 63 seqfil = seqfile 61 endif 64 endif 62 65 call redseq 63 64 write (*,*) 'File with sequence is ', seqfil(1:iendst(seqfil)) 65 66 67 write (logString, *) 'File with sequence is ', 68 & seqfil(1:iendst(seqfil)) 69 66 70 ! --------------------------------- read & assemble data from libraries 67 71 ! initial coordinates, interaction lists 68 72 69 73 ntl = ntlml 70 74 do i=ontlml, ntl 71 75 72 76 call getmol(i) ! assemble data from libraries 73 77 74 78 do j=1,6 ! initialize global parameters 75 79 gbpr(j,i)=0.d0 76 80 enddo 77 81 78 82 call bldmol(i) ! co-ordinates 79 83 80 84 ntlml = i 81 85 call addend(i,grpn,grpc) ! modify ends 82 86 call setmvs(i) ! determine sets of moving atoms for given variables 83 87 call mklist(i) ! compile lists of interaction partners 84 88 85 89 enddo 86 90 87 91 ! --------------------------- Read the initial conformation if necessary 88 if(readFromStdin) then 89 write ( *,'(a,$)') ' file with VARIABLES:'90 ! 92 if(readFromStdin) then 93 write (logString, '(a,$)') ' file with VARIABLES:' 94 ! 91 95 varfil=' ' 92 96 read(*,'(a)',end=2,err=2) varfil … … 96 100 l=iendst(varfil) 97 101 if (l.gt.0.and.varfil.ne.' ') then 98 write ( *,'(1x,a,/)') varfil(1:l)102 write (logString, '(1x,a,/)') varfil(1:l) 99 103 lunvar=13 100 104 101 105 call redvar ! get vars. and rebuild 102 103 endif 104 105 2 write(*,*) ' ' 106 107 ireg = 0 108 109 else ! =========================== from PDB 110 if (iendst(seqfile).le.1) then 111 3 write (*,'(/,a,$)') ' PDB-file:' 112 seqfil=' ' 113 read (*,'(a)',err=3) seqfil 114 else 115 seqfil = seqfile 116 endif 117 write (*,*) 'PDB structure ',seqfil(1:iendst(seqfil)) 118 print *, 'calling readpdb with ',seqfile 119 call pdbread(seqfil,ier) 120 121 if (ier.ne.0) stop 122 123 call pdbvars() 124 125 ireg = 1 126 127 endif 128 129 ! If Lund force field is in use, keep omega angles fixed 130 if (ientyp.eq.2) then 131 do iv=1,nvrml(ntlml) 132 if ((nmvr(iv)(1:2).eq.'om')) then 133 vlvr(iv)=pi 134 toat(iatvr(iv))=pi 135 fxvr(iv)=.true. 136 print *, 'Fixed variable ',iv,nmvr(iv),vlvr(iv) 137 endif 138 enddo 139 endif 106 107 endif 108 109 2 write (logString, *) ' ' 140 110 141 111 ! -------------------- get: nvr,idvr, vlvr, olvlvr … … 161 131 enddo 162 132 133 ireg = 0 134 135 else ! =========================== from PDB 136 if (iendst(seqfile).le.1) then 137 3 write (logString, '(/,a,$)') ' PDB-file:' 138 seqfil=' ' 139 read (*,'(a)',err=3) seqfil 140 else 141 seqfil = seqfile 142 endif 143 write (logString, *) 'PDB structure ',seqfil(1:iendst(seqfil)) 144 print *, 'calling readpdb with ',seqfile 145 call pdbread(seqfil,ier) 146 147 if (ier.ne.0) stop 148 149 call pdbvars() 150 151 ireg = 1 152 153 endif 154 163 155 ! -------------------------- set var. amplitudes for simulations 164 156 165 157 do i=1,ivrml1(ntlml)+nvrml(ntlml)-1 166 158 167 159 if (ityvr(i).eq.3.and..not.fxvr(i)) then ! torsion 168 160 169 161 navr = nmvr(i) 170 162 171 163 ir = nursvr(i) 172 164 nars = seq(ir) 173 165 174 166 if ( navr(1:2).eq.'om' 175 167 176 168 & .or.nars(1:3).eq.'arg'.and.(navr(1:2).eq.'x5' 177 169 & .or.navr(1:2).eq.'x6') 178 170 179 171 & .or.(nars(1:3).eq.'asn'.or.nars(1:3).eq.'asp') 180 172 & .and.navr(1:2).eq.'x3' 181 173 182 174 & .or.(nars(1:3).eq.'gln'.or.nars(1:3).eq.'glu') 183 175 & .and.navr(1:2).eq.'x4' 184 176 185 177 & ) then 186 178 187 179 ! axvr(i) = pi/9.d0 ! 20 deg. 188 180 axvr(i) = pi2 ! Trying out 360 deg. for these as well 189 181 190 182 else 191 183 axvr(i) = pi2 ! 360 deg. 192 184 endif 193 185 194 186 else 195 187 axvr(i) = 0.d0 196 188 endif 197 189 198 190 enddo ! vars. 199 191 200 192 ! --------------------- initialize solvation pars. if necessary 201 193 202 194 if (itysol.ne.0) then 203 195 204 196 i1=iatrs1(irsml1(1)) ! 1st atom of 1st molecule 205 197 i2=iatrs2(irsml2(ntlml)) ! last atom of last molecule 206 198 207 199 its = iabs(itysol) 208 200 209 201 do i=i1,i2 ! all atoms 210 202 it=ityat(i) 211 203 sigma(i)=coef_sl(its,it) 212 204 rvdw(i) =rad_vdw(its,it) 213 205 214 206 if (nmat(i)(1:1).ne.'h') rvdw(i)=rvdw(i)+rwater 215 207 216 208 enddo 217 209 218 210 endif 219 211 ! Initialize calpha array 220 212 do i=ontlml, ntlml 221 call c_alfa(i,1) 213 call c_alfa(i,1) 222 214 enddo 223 215 224 216 ! Initialize arrays used in the BGS update 225 call init_lund() 217 call init_lund() 226 218 return 227 219 end 228 229 220 221 -
main.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 16 16 ! $Id: main.f 334 2007-08-07 09:23:59Z meinke $ 17 17 ! ************************************************************** 18 18 19 19 program main 20 20 21 21 include 'INCL.H' 22 22 include 'INCP.H' 23 double precision eps, temp, e, energy 24 25 integer maxloglevel, logfileunit, iabin, imin, maxit, nequi 26 integer nsweep, nmes, ncalls, nacalls 27 23 28 common/updstats/ncalls(5),nacalls(5) 24 29 character*80 libdir, seqfile, varfile … … 32 37 ! Directory for SMMP libraries 33 38 ! Change the following directory path to where you want to put SMMP 34 ! libraries of residues. 39 ! libraries of residues. 35 40 libdir='./SMMP/' 36 41 37 42 ! Set the maximum log level. The larger the number the more detailed 38 43 ! the log. … … 48 53 ientyp = 0 49 54 ! 0 => ECEPP2 or ECEPP3 depending on the value of sh2 50 ! 1 => FLEX 55 ! 1 => FLEX 51 56 ! 2 => Lund force field 52 57 ! 3 => ECEPP with Abagyan corrections … … 70 75 iabin = 1 ! =0: read from PDB-file 71 76 ! =1: ab Initio from sequence (& variables) 72 seqfile='EXAMPLES/enkefa.seq' 73 varfile='EXAMPLES/enkefa.ann' 74 ! varfile = ' ' 75 77 seqfile='polyq.seq' 78 ! seqfile='polyA.pdb' 79 varfile='polyq.var' 80 varfile = ' ' 81 76 82 ntlml = 0 77 83 write (*,*) 'Solvent: ', itysol 78 84 ! Initialize random number generator. 79 85 call sgrnd(31433) 80 86 81 87 if (itysol.eq.0.and.ientyp.eq.3) then 82 88 print *,'Can not use Abagyan entropic corrections without ' … … 87 93 call init_molecule(iabin,grpn,grpc,seqfile,varfile) 88 94 89 ! Decide if and when to use BGS, and initialize Lund data structures 95 ! Decide if and when to use BGS, and initialize Lund data structures 90 96 bgsprob=0.75 ! Prob for BGS, given that it is possible 91 ! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 92 ! 2 => temperature dependent choice 97 ! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 98 ! 2 => temperature dependent choice 93 99 upchswitch=1 94 100 rndord=.true. … … 96 102 if (ientyp.eq.2) call init_lundff 97 103 if (ientyp.eq.3) call init_abgn 98 104 99 105 100 106 ! ======================================== Add your task down here … … 103 109 maxit = 15000 ! maximum number of iterations in minimization 104 110 eps = 1.0d-7 ! requested precision 105 call minim(imin, maxit, eps)106 call outvar(0, ' ')111 ! call minim(imin, maxit, eps) 112 ! call outvar(0, ' ') 107 113 ! To do a canonical Monte Carlo simulation uncomment the lines below 108 ! nequi = 100 109 ! nsweep = 50000 110 ! nmes = 10 111 ! temp = 300.0 112 ! lrand = .true. 114 nequi = 100 115 nsweep = 50000 116 nmes = 10 117 temp = 300.0 118 lrand = .true. 119 E = energy() 120 write(*,*) E, eyel,eyvw,eyhb,eyvr 121 call outpdb(1, "polyA.pdb") 113 122 ! Canonical Monte Carlo 114 ! call canon(nequi, nsweep, nmes, temp, lrand)123 ! call canon(nequi, nsweep, nmes, temp, lrand) 115 124 116 125 ! For simulated annealing uncomment the lines below … … 118 127 ! tmax = 500.0 119 128 ! call anneal(nequi, nsweep, nmes, tmax, tmin, lrand); 120 ! ======================================== End of main 129 ! ======================================== End of main 121 130 end -
main_p.f
r6650a56 r32289cd 1 1 ! ************************************************************** 2 ! 2 ! 3 3 ! This file contains the main (PARALLEL TEMPERING JOBS ONLY, 4 4 ! FOR SINGULAR PROCESSOR JOBS USE main) 5 ! 5 ! 6 6 ! This file contains also the subroutine: p_init_molecule 7 ! 7 ! 8 8 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 9 ! Shura Hayryan, Chin-Ku 9 ! Shura Hayryan, Chin-Ku 10 10 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 11 11 ! Jan H. Meinke, Sandipan Mohanty 12 ! 12 ! 13 13 ! CALLS init_energy,p_init_molecule,partem_p 14 ! 14 ! 15 15 ! ************************************************************** 16 16 program pmain … … 21 21 include 'mpif.h' 22 22 23 double precision startwtime, group_world, error, endwtime 24 25 integer ierr, num_proc, iabin, nequi, nswp, nmes, nsave, ifrm, j 26 integer i, nml, nresi, my_pt_rank, ncalls, nacalls 27 23 28 character*80 libdir 24 29 character*80 in_fil,ou_fil,filebase, varfile … … 28 33 logical newsta 29 34 30 !c Number of replicas 35 !c Number of replicas 31 36 integer num_replica 32 37 !c Number of processors per replica … … 62 67 63 68 ! =================================================== Energy setup 64 libdir='SMMP/' 69 libdir='SMMP/' 65 70 ! Directory for SMMP libraries 66 71 … … 71 76 ientyp = 0 72 77 ! 0 => ECEPP2 or ECEPP3 depending on the value of sh2 73 ! 1 => FLEX 78 ! 1 => FLEX 74 79 ! 2 => Lund force field 75 80 ! 3 => ECEPP with Abagyan corrections 76 ! 81 ! 77 82 78 83 sh2=.false. ! .true. for ECEPP/2; .false. for ECEPP3 … … 106 111 ! temperatures must have at least as many 107 112 ! entries 108 nequi=10 ! Number of MC sweeps before measurements 109 ! and replica exchanges are started 113 nequi=10 ! Number of MC sweeps before measurements 114 ! and replica exchanges are started 110 115 nswp=500000 ! Number of sweeps 111 116 nmes=10 ! Interval for measurements and replica exchange 112 117 nsave=1000 ! Not used at the moment 113 114 switch = -1 ! How should the configuration be 118 119 switch = -1 ! How should the configuration be 115 120 ! initialized? 116 ! -1 stretched chain 121 ! -1 stretched chain 117 122 ! 0 don't do anything 118 123 ! 1 initialize each angle to a random value 119 124 120 125 ifrm=0 121 126 ntlml = 0 122 127 123 ! Decide if and when to use BGS, and initialize Lund data structures 128 ! Decide if and when to use BGS, and initialize Lund data structures 124 129 bgsprob=0.6 ! Prob for BGS, given that it is possible 125 ! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 126 ! 2 => temperature dependent choice 130 ! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 131 ! 2 => temperature dependent choice 127 132 upchswitch=1 128 133 rndord=.true. 129 134 ! ================================================================= 130 135 ! Distribute nodes to parallel tempering tasks 131 ! I assume that the number of nodes available is an integer 136 ! I assume that the number of nodes available is an integer 132 137 ! multiple n of the number of replicas. Each replica then gets n 133 138 ! processors to do its energy calculation. … … 139 144 ! could just use i * num_ppr but this way it's more flexible. 140 145 j = 0 141 do i = 1, num_replica 142 ranks(i) = j 146 do i = 1, num_replica 147 ranks(i) = j 143 148 proc_range(1) = j 144 149 proc_range(2) = j + num_ppr - 1 … … 146 151 call mpi_group_range_incl(group_world, 1, proc_range, group(i) 147 152 & ,error) 148 write (*,*) "Assigning rank ", j, proc_range, 153 write (*,*) "Assigning rank ", j, proc_range, 149 154 & "to group", group(i) 150 155 call flush(6) … … 167 172 call mpi_group_incl(group_world, num_replica, ranks, group_partem, 168 173 & error) 169 call mpi_comm_create(mpi_comm_world, group_partem, partem_comm, 174 call mpi_comm_create(mpi_comm_world, group_partem, partem_comm, 170 175 & error) 171 176 … … 184 189 varfile = 'EXAMPLES/1bdd.var' 185 190 call init_molecule(iabin, grpn, grpc,in_fil,varfile) 186 else 191 else 187 192 filebase = "conf_0000.var" 188 193 call init_molecule(iabin, grpn, grpc,in_fil, … … 220 225 nml = 1 221 226 call distributeWorkLoad(no, nml) 222 227 223 228 call partem_p(num_replica, nequi, nswp, nmes, nsave, newsta, 224 229 & switch, rep_id, partem_comm) -
metropolis.f
r6650a56 r32289cd 20 20 include 'INCP.H' 21 21 include 'incl_lund.h' 22 23 double precision eol, energy, grnd, dv, addang, enw, delta, dummy 24 double precision rd, ex, acz, eol1 25 26 integer nsw, iupstate, iupt, ivar, jv, iml, i, ncalls, nacalls 27 22 28 external dummy 23 29 ! common/bet/beta … … 178 184 179 185 subroutine accanalyze(iuptype,iupdstate) 186 integer ncalls, nacalls, iupdstate, iuptype 180 187 common/updstats/ncalls(5),nacalls(5) 181 188 ncalls(5)=ncalls(5)+1 … … 202 209 integer iiii 203 210 double precision bbbb,curprob, grnd,up2bmax,up2bmin 211 204 212 common/updtparam/up2bmax,up2bmin 205 213 curprob=(bbbb-up2bmin)/(up2bmax-up2bmin) … … 211 219 end 212 220 block data updtchs 221 integer ncalls, nacalls 213 222 double precision up2bmax, up2bmin 214 223 common/updtparam/up2bmax,up2bmin -
mincjg.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 40 40 ! .................................................................... 41 41 42 implicit real*8 (a-h,o-z) 43 implicit integer*4 (i-n) 42 double precision amf, tol, eps, gsqrd, fmin, f, t, g, ga, d, xa, x 43 double precision gsq2, gnew, dfpr, stmin, fuit, gdit, gt, gdmi 44 double precision sbound, stepch, step, wo, sum, fch, acur, ddspln 45 double precision gspln, beta, gama, yt, dt, gamden, di 46 47 integer mxfcon, maxlin, mxn, ier, iter, iterfm, iterrs, nfun 48 integer nfopt, i, n, nfbeg, iretry, maxfun 44 49 45 50 parameter (AMF = 10.d0, … … 102 107 103 108 2 step = stmin + stepch 104 wo = 0.d0 109 wo = 0.d0 105 110 106 111 do i=1,n … … 182 187 183 188 if ( fch .ne. 0.d0 ) ddspln = ddspln + (wo + wo) / stepch 184 189 185 190 if ( gdmi .eq. 0.d0 ) goto 6 186 191 -
nursvr.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 18 18 ! ........................................................... 19 19 include 'INCL.H' 20 21 integer i, ifirs, ivr, j 20 22 21 23 do i=ntlml,1,-1 … … 45 47 include 'INCL.H' 46 48 49 integer i, ifirs, ilars, iat, j 50 47 51 do i=1,ntlml 48 52 -
opereg.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 33 33 include 'INCL.H' 34 34 include 'INCP.H' 35 36 double precision xfvr, yfvr, zfvr, xfrvr, yfrvr, zfrvr, x1, y1, z1 37 double precision a, sa, ca, xk, yk, zk, xi, yi, zi, xji, yji, zji 38 double precision dx, dy, dz, x, y, z, xfat, yfat, zfat, xfrat 39 double precision yfrat, zfrat, xb, yb, zb, ex, ey, ez, xfiv, yfiv 40 double precision zfiv, xfriv, yfriv, zfriv 41 42 integer ntlvr, nml, ix2, ifivr, ilavr, i, ii, j, i1s, i1a, io, iv 43 integer it, ia, ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, i14 35 44 36 45 dimension xfat(mxat),yfat(mxat),zfat(mxat), … … 100 109 101 110 eyrg = eyrg + xji**2 + yji**2 + zji**2 ! The regularization energy is just 102 ! the sum over the atom distances 111 ! the sum over the atom distances 103 112 ! squared. 104 113 … … 303 312 include 'INCL.H' 304 313 314 double precision del, pro, gdn, enyreg 315 316 integer ii, nml, i 317 305 318 parameter (del=1.d-7) 306 319 … … 339 352 340 353 include 'INCL.H' 354 355 double precision del, vlvrx, ovr, eynw, enyreg, gdn, gda 356 357 integer i, it, iv, nml 341 358 342 359 parameter (del=1.d-6) -
opesol.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 39 39 integer ntlvr, nml, ix2, ifivr, ilavr, i, i1s, i1a, io, iv, it, ia 40 40 integer ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, j, i14 41 41 42 42 dimension xfat(mxat),yfat(mxat),zfat(mxat), 43 43 & xfrat(mxat),yfrat(mxat),zfrat(mxat), … … 251 251 252 252 include 'INCL.H' 253 253 254 254 double precision del, vlvrx, ovr, eynw, esolan, gda, gdn 255 255 -
outpdb.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 24 24 25 25 include 'INCL.H' 26 27 double precision occ, bva 28 29 integer i0, i9, nml, im1, im2, ibegst, iout, iat, iml, irs, ifirs 30 integer ifiat, nrs, i, j, iendst, nfi, ibd, jj, nbd 26 31 27 32 dimension ibd(4) … … 99 104 call toupst(atnm) 100 105 call toupst(res) 101 106 102 107 linout = ' ' 103 108 write (linout,1) linty,iat,atnm,res(1:3),chid,irs,cdin, … … 170 175 write (linout,'(a3)') 'END' 171 176 write(iout,'(a80)') linout 172 177 173 178 close(iout) 174 179 -
outvar.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty 9 9 ! 10 10 ! ************************************************************** 11 11 12 12 subroutine outvar(nml,fileName) 13 13 14 14 !-------------------------------------------------------------------- 15 ! Output of variables of the current protein conformation 16 !- 15 ! Output of variables of the current protein conformation 16 !- 17 17 ! nml != 0 : molecule index 18 18 ! nml = 0 : all molecules … … 24 24 25 25 include 'INCL.H' 26 27 integer nml, ibegst, iout 26 28 27 29 character*(*) fileName … … 49 51 50 52 include 'INCL.H' 53 54 integer nml, im1, im2, iml, i, iout, ibegst, ifivr, is, nursvr, it 51 55 52 56 character mlfd*7,strg(6)*17 -
partem_p.f
r6650a56 r32289cd 43 43 44 44 integer ifrrm, nmes, nswp, num_rep, i, j, nresi, iold, inode 45 integer intem, iv, jold, idum1, idum2, idum3 , mpi_integer46 integer mpi_comm_world, ierr, mpi_double_precision, nsw, nequi45 integer intem, iv, jold, idum1, idum2, idum3 46 integer ierr, nsw, nequi 47 47 integer nml, nhel, mhel, nbet, mbet, mhb, imhb, nctot, ncnat 48 integer mpi_comm_null,k1, k, nu, no1, in, jn48 integer k1, k, nu, no1, in, jn 49 49 50 50 dimension eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC), … … 225 225 call flush(6) 226 226 endif 227 if(mod(iold,nmes).eq.0 ) then227 if(mod(iold,nmes).eq.0.and.myrank.eq.0) then 228 228 if ((rep_id + 1).eq.trackID.and.myrank.eq.0) then 229 229 frame = iold /nmes -
partem_s.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 10 10 ! 11 11 ! ************************************************************** 12 12 13 13 subroutine partem_s(num_rep, nequi, nswp, nmes, nsave, newsta, 14 14 & switch) … … 20 20 ! 21 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 22 29 external can_weight 23 30 ! TODO Store global coordinates in pgbpr … … 44 51 ! ipoi: Points to replica ipoi(k) which is currently 45 52 ! at inverse temperature pbe(k) 46 ! eol: energy of each replica 53 ! eol: energy of each replica 47 54 ! acc: accepatance rate of each replica 48 55 ! … … 55 62 allocate(temp(num_rep), pbe(num_rep), eol(num_rep)) 56 63 allocate(acc(num_rep), ipoi(num_rep)) 57 64 58 65 if (.not.(allocated(coor_G).and.allocated(temp) 59 66 & .and. allocated(pbe).and. allocated(eol) … … 77 84 ipoi(i) = 0 78 85 end do 79 86 80 87 81 88 ! READ IN TEMPERATURES … … 90 97 91 98 open(13,file='time.d',status='unknown') 92 99 93 100 if(.not.newsta) then 94 101 ! READ START Values … … 122 129 pgbpr(j, nml, i) = gbpr(j, nml) 123 130 end do 124 end do 125 end do 126 131 end do 132 end do 133 127 134 ! Equilibrization for each replica (No replica exchange move) 128 135 do k=1,num_rep … … 141 148 CALL METROPOLIS(energ,acz,can_weight) 142 149 end do 143 write(*,*) 'Start energy after equilibration for replica:', 150 write(*,*) 'Start energy after equilibration for replica:', 144 151 & k, energ 145 152 do i=1,nvr … … 155 162 end do 156 163 end if 157 164 158 165 ! Now begins the simulation with Multiple Markov Chains 159 166 iswitch = 1 … … 176 183 end do 177 184 CALL METROPOLIS(energ,acz,can_weight) 178 ! 185 ! 179 186 if(mod(nsw,nmes).eq.0) then 180 187 ! Measure and store here all quantities you want to analyse later … … 195 202 end do 196 203 end do 197 204 198 205 eol(k) = energ 199 206 acc(k) = acz … … 214 221 j=i+1 215 222 if(i.eq.num_rep) j=1 216 in=ipoi(i) 223 in=ipoi(i) 217 224 jn=ipoi(j) 218 225 delta=-pbe(i)*eol(jn)-pbe(j)*eol(in) -
pdbread.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 21 21 ! ...................................................... 22 22 23 implicit real*8 (a-h,o-z)24 implicit integer*4 (i-n)25 26 23 include 'INCP.H' 24 25 double precision cor 26 27 integer ier, l, iendst, lunpdb, io, iopfil, iat, i 27 28 28 29 ! -------------------------- input … … 47 48 lunpdb = 99 48 49 else 49 write (*,'(a)') 50 write (*,'(a)') 50 51 & ' pdbread> empty file name to read pdb-structure' 51 52 … … 56 57 57 58 if (io.le.0) then 58 write (*,'(a,/,a)') 59 write (*,'(a,/,a)') 59 60 & ' pdbread> ERROR opening file to read pdb-structure: ', 60 61 & pdbfil(1:iendst(pdbfil)) … … 69 70 70 71 if ( line(17:17).ne.' ' ) then 71 write (*,'(a,/,a,/,a,/,2a)') 72 write (*,'(a,/,a,/,a,/,2a)') 72 73 & ' pdbread> found alternate atom location: ', 73 74 & ' !', … … 85 86 86 87 if ((natp+1).gt.MXATP) then 87 write (*,'(a,i5,a,/,a)') 88 write (*,'(a,i5,a,/,a)') 88 89 & ' pdbread> >MXATP (',MXATP,') ATOM lines in PDB file ', 89 90 & pdbfil(1:iendst(pdbfil)) … … 96 97 97 98 if ((nchp+1).gt.MXCHP) then 98 write (*,'(a,i3,a,/,a)') 99 write (*,'(a,i3,a,/,a)') 99 100 & ' pdbread> >MXCHP (',MXCHP,') chains in PDB file ', 100 101 & pdbfil(1:iendst(pdbfil)) … … 105 106 106 107 if ((nrsp+1).gt.MXRSP) then 107 write (*,'(a,i3,a,/,a)') 108 write (*,'(a,i3,a,/,a)') 108 109 & ' pdbread> >MXRSP (',MXRSP,') residues in PDB file ', 109 110 & pdbfil(1:iendst(pdbfil)) … … 140 141 141 142 if ((nrsp+1).gt.MXRSP) then 142 write (*,'(a,i3,a,/,a)') 143 write (*,'(a,i3,a,/,a)') 143 144 & ' pdbread> >MXRSP (',MXRSP,') residues in PDB file ', 144 145 & pdbfil(1:iendst(pdbfil)) … … 171 172 goto 1 172 173 173 2 write (*,'(a,/,a,/,2a)') 174 2 write (*,'(a,/,a,/,2a)') 174 175 & ' pdbread> ERROR reading ATOM line ', 175 176 & line(:l), … … 194 195 else 195 196 196 write (*,'(a,/,a)') 197 write (*,'(a,/,a)') 197 198 & ' pdbread> NO atom coordinates selected from file ', 198 199 & pdbfil(1:iendst(pdbfil)) … … 226 227 include 'INCL.H' 227 228 include 'INCP.H' 229 230 double precision dihedr, rm, av1, av2, rmsd, h, x, y, z 231 232 integer nml, nrs, nc, irb, ire, irs, i, ii, it, i1, i2, i3, i4, j1 233 integer j2, j3, j4, inew, j, ix 228 234 229 235 character res*3 … … 451 457 include 'INCP.H' 452 458 459 integer irs, nml, i1, i2, iat, ix, i 460 453 461 character*4 atm 454 462 … … 476 484 endif 477 485 enddo 478 486 479 487 ! write(*,'(8a)') ' pdbvars> ',atm,' not found in ' 480 488 ! # ,chnp(nc),' ',rsidp(irs),' ',rsnmp(irs) … … 484 492 1 ixatp(iat)=ix 485 493 486 enddo ! SMMP atoms of 'irs' 494 enddo ! SMMP atoms of 'irs' 487 495 enddo ! residues 488 496 … … 494 502 include 'INCL.H' 495 503 504 double precision tol, h1, h2, h3, x1, x2, x3, d, z1, z2, z3, y1 505 double precision y2, y3, yp1, yp2 506 507 integer i1, nml, i2, i3, i 508 496 509 parameter (TOL = 1.d-12) 497 510 … … 499 512 ! -> determine global parameters: shifts dX,dY,dZ 500 513 ! & angles alpha,beta,gamma [rad], put into 'gbpr' 501 ! 514 ! 502 515 ! CALLS: none 503 516 ! … … 547 560 y2=z3*x1-z1*x3 548 561 y3=z1*x2-z2*x1 549 562 550 563 if ( ( 1.d0 - abs(y3) ) .gt. TOL ) then ! ============ |beta| < PI/2 551 564 -
regul.f
r6650a56 r32289cd 4 4 ! 5 5 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 6 ! Shura Hayryan, Chin-Ku 6 ! Shura Hayryan, Chin-Ku 7 7 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 8 8 ! Jan H. Meinke, Sandipan Mohanty … … 10 10 ! ************************************************************** 11 11 12 12 13 13 subroutine regul(nml, iter, nsteps, acc) 14 14 … … 26 26 include 'INCL.H' 27 27 include 'INCP.H' 28 28 29 29 !f2py intent(in) nml 30 30 !f2py intent(in) iter 31 31 !f2py intent(in) nsteps 32 32 !f2py intent(in) acc 33 34 double precision acc, rm, av1, av2, rmsd, dn 35 36 integer nsteps, nml, nrs, i, n, iter, it 33 37 34 38 dimension rm(3,3),av1(3),av2(3) … … 60 64 61 65 62 do i = ivrml1(nml),nvrml(nml) 66 do i = ivrml1(nml),nvrml(nml) 63 67 fxvro(i) = fxvr(i) ! save 64 68 if (isrfvr(i)) fxvr(i) = .true. ! fix vars. defined in ref.str. … … 77 81 call cnteny(nml) 78 82 79 do i = ivrml1(nml),nvrml(nml) 83 do i = ivrml1(nml),nvrml(nml) 80 84 fxvr(i) = fxvro(i) ! restore 81 85 enddo ! vars. … … 86 90 wtey = 0.d0 87 91 88 n=iter 92 n=iter 89 93 dn=1.d0/dble(n) 90 94 -
rgyr.f
r6650a56 r32289cd 25 25 ! 26 26 include 'INCL.H' 27 28 double precision dn, dnp, dnh, dx, dxp, dy, dyp, dyh, dz, dzp, dzh 29 double precision d2, d2p, d2h, xi, yi, zi, dxh, rg2, rg2p, rg2h 30 double precision rgy, ee 31 32 integer nml, nml1, nml2, nat, i, i1, i2 33 27 34 !f2py intent(in) nml 28 35 !f2py intent(out) rgy -
rmsdfun.f
r6650a56 r32289cd 5 5 ! 6 6 ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann, 7 ! Shura Hayryan, Chin-Ku 7 ! Shura Hayryan, Chin-Ku 8 8 ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann, 9 9 ! Jan H. Meinke, Sandipan Mohanty … … 26 26 ! 27 27 ! Input 28 double precision xrf, yrf, zrf, rm, av1, av2, rssd 29 30 integer nml, ir1, ir2, ixat, isl 31 28 32 dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp) 29 33 ! Local … … 57 61 ! 58 62 ! NB uncomment last lines in 'fitmol' to return coordinates 59 ! in 'x2' after fitting the ref. str. onto SMMP structure 63 ! in 'x2' after fitting the ref. str. onto SMMP structure 60 64 ! ---------------------------------------------------------------- 61 65 62 66 include 'INCL.H' 63 67 include 'INCP.H' 68 69 double precision x1, x2, xrf, yrf, zrf, rm, av1, av2, rmsd 70 71 integer nml, nr, na, n, im, ir, ia, ir1, ir2, ix, ixat, isl 72 64 73 65 74 !-------------------------------------------------------- input … … 147 156 ! ....................................................... 148 157 !f2py intent(out) rmsd 149 158 150 159 include 'INCL.H' 151 160 ! implicit real*8 (a-h,o-z) 152 161 ! implicit integer*4 (i-n) 153 162 154 163 ! ------------------------------------------- input/output 164 double precision dn, a1, a2, x1, x2, q, dm, dp, dxm, dym, dzm, dxp 165 double precision dyp, dzp, e, v, em, rmsd, rm 166 167 integer n, i, j, ndim4, im 168 155 169 dimension x1(3,mxat),x2(3,mxat) 156 170 ! -------------------------------------------------- local … … 270 284 !f2py intent(out) d 271 285 !f2py intent(out) v 286 integer nmax 287 272 288 parameter (NMAX=500) 273 289 274 290 integer n,nrot,i,ip,iq,j 275 291 … … 398 414 ! 399 415 !------------------------------------------------------------------------------ 400 ! Reads in pdb-file 'string' into INCP.H and initalizes 416 ! Reads in pdb-file 'string' into INCP.H and initalizes 401 417 ! the files that 'rmdsopt' needs to calculate the rmsd 402 418 ! of a configuration with the pdb-configuration … … 409 425 include 'INCP.H' 410 426 427 integer i, nml, ier 428 411 429 character string*(*) 412 430 413 431 if(string.eq.'smmp') then 414 432 ! … … 428 446 ! 429 447 call pdbread(string,ier) 430 if(ier.ne.0) stop 448 if(ier.ne.0) stop 431 449 call atixpdb(nml) 432 450 ! -
twister.f
r6650a56 r32289cd 47 47 * May 03, 2000 48 48 * 49 * compile with: f90 -i32 -dp -o mt19937 mt19937.f 49 * compile with: f90 -i32 -dp -o mt19937 mt19937.f 50 50 * Changes: 51 51 * Do not use UMASK as parameter due to restrictions of CF90 … … 55 55 * 56 56 * Now the T3E produces the same random number as a RS6000 running AIX 57 * 57 * 58 58 * from other programs sgrnd() must be called with integer*4 59 59 * … … 64 64 subroutine sgrnd(seed) 65 65 * 66 implicit integer(a-z) 66 integer seed, tmaskb, tmaskc, tshftu, tshfts, tshftt 67 integer tshftl, y 68 69 integer N, N1, mti, mt, m, mata, lmask, mag01 70 integer kk 67 71 * 68 72 * Period parameters … … 88 92 89 93 block data twbloks 94 integer n,n1,mt,mti 90 95 parameter(N = 624) 91 96 parameter(N1 = N+1) … … 99 104 double precision function grnd() 100 105 * 101 implicit integer(a-z) 106 integer tmaskb, tmaskc, tshftu, tshfts, tshftt, tshftl, y 107 108 integer n, n1, mti, mt, m, mata, lmask, mag01, kk 109 102 110 * 103 111 * Period parameters -
utilities.f
r6650a56 r32289cd 8 8 9 9 !! Calculate the best way to distribute the work load across processors. 10 ! It calculates the average number of interactions and then tries to 10 ! It calculates the average number of interactions and then tries to 11 11 ! assign a number of interactions to each processor that is as close 12 12 ! as possible to the average. The routine should be called once for … … 21 21 include 'INCL.H' 22 22 23 integer i1ms, io, iv, i2ms, ms 24 23 25 integer num_ppr, nml 24 26 integer idxOfFirstVariable, idxOfLastVariable … … 26 28 integer totalct, irank, itarget 27 29 double precision ipps 28 30 29 31 if (nml.eq.0) then 30 32 idxOfFirstVariable = ivrml1(1) … … 36 38 end do 37 39 end do 38 else 40 else 39 41 idxOfFirstVariable = ivrml1(nml) 40 42 idxOfLastVariable = ivrml1(nml) + nvrml(nml) - 1 … … 44 46 end do 45 47 end if 46 48 47 49 isum = 0 48 50 do io = idxOfLastVariable, idxOfFirstVariable, - 1 … … 52 54 do ms = i1ms, i2ms 53 55 do at = latms1(ms), latms2(ms) 54 do ivw=ivwat1(at),ivwat2(at) 56 do ivw=ivwat1(at),ivwat2(at) 55 57 do j=lvwat1(ivw),lvwat2(ivw) 56 58 isum = isum + 1 57 59 end do 58 60 end do 59 do i14=i14at1(at),i14at2(at) 61 do i14=i14at1(at),i14at2(at) 60 62 isum = isum + 1 61 63 end do … … 66 68 write (*,*) "Total number of interactions:", isum 67 69 write (*,*) "Average # of interactions per processor", ipps 68 70 69 71 totalct = 0 70 72 irank = 1 … … 72 74 if (nml.eq.0) then 73 75 i1ms = imsml1(ntlml)+ nmsml(ntlml) 74 else 76 else 75 77 i1ms = imsml1(nml)+ nmsml(nml) 76 78 end if 77 79 do io = idxOfLastVariable, idxOfFirstVariable, - 1 78 isum = 0 80 isum = 0 79 81 iv = iorvr(io) 80 82 i2ms = i1ms - 1 … … 83 85 do at = latms1(ms), latms2(ms) 84 86 atct = atct + 1 85 do ivw=ivwat1(at),ivwat2(at) 87 do ivw=ivwat1(at),ivwat2(at) 86 88 do j=lvwat1(ivw),lvwat2(ivw) 87 89 isum = isum + 1 88 90 end do 89 91 end do 90 do i14=i14at1(at),i14at2(at) 92 do i14=i14at1(at),i14at2(at) 91 93 isum = isum + 1 92 94 end do … … 113 115 114 116 end subroutine distributeWorkLoad 115 117 116 118 !----------------------------------------------------------------------- 117 119 ! The function fileNameMP takes a template of a file name in the … … 128 130 ! \endcode 129 131 ! will output base_0011.dat. 130 ! 132 ! 131 133 ! @param base the base file name, e.g., base_0000.dat. 132 134 ! @param i1 index of the first character that may be replaced 133 135 ! @param i2 index of the last character that may be replaced 134 136 ! @param rank the number that should be inserted into the file name. 135 ! 137 ! 136 138 ! @return file name for rank 137 139 !----------------------------------------------------------------------- … … 165 167 write(fileNameMP(i2-5:i2), '(i6)') rank 166 168 endif 167 end function fileNameMP 169 end function fileNameMP 168 170 ! End fileNameMP 169 171 … … 172 174 ! Add messages to log. This routine takes the log (debugging) mes- 173 175 ! sages and writes them to the log file if the log level is less or 174 ! equal to the maximum log level given by the global variable 176 ! equal to the maximum log level given by the global variable 175 177 ! MAXLOGLEVEL. 176 178 ! 177 179 ! @author Jan H. Meinke 178 180 ! 179 ! @param loglevel level at which this message should be added to 181 ! @param loglevel level at which this message should be added to 180 182 ! the log. 181 183 ! @param message message to be written to the log. 182 ! @param rank global rank of this node if running an MPI job zero 184 ! @param rank global rank of this node if running an MPI job zero 183 185 ! otherwise. 184 186 !---------------------------------------------------------------------- 185 187 subroutine addLogMessage(loglevel, message, rank) 186 188 189 integer maxloglevel, logfileunit 190 187 191 integer :: loglevel, rank 188 192 character(LEN=*) :: message 189 193 190 194 if (loglevel <= MAXLOGLEVEL) then 191 195 write(LOGFILEUNIT, *) message 192 196 end if 193 197 194 198 end subroutine addLogMessage -
zimmer.f
r6650a56 r32289cd 40 40 !f2py intent(in) nresi 41 41 character*1 zim 42 double precision xphi, xpsi 43 44 integer j, nresi, i, iv, nres, nursvr 42 45 43 46
Note:
See TracChangeset
for help on using the changeset viewer.