From f061b417840ab72282fafafd36d90a33e7de34c9 Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Fri, 23 May 2014 15:41:33 -0400 Subject: [PATCH] ported isokinetic and stochastic isokinetic dynamics from header-tinker to module-tinker --- source/analysis.f | 52 ++--- source/angang.f | 1 + source/angbnd.f | 1 + source/angles.f | 4 +- source/attach.f | 12 +- source/bar.f | 4 +- source/bath.f | 13 ++ source/bitors.f | 3 - source/bonds.f | 4 +- source/cluster.f | 2 +- source/connolly.f | 349 +++++++++++----------------- source/diffuse.f | 55 ++--- source/domega.f | 53 ++--- source/dynamic.f | 14 +- source/eangang2.f | 4 - source/embed.f | 12 +- source/eopbend2.f | 4 - source/esolv1.f | 6 +- source/faces.f | 112 +++++---- source/final.f | 125 +--------- source/flatten.f | 6 - source/getkey.f | 1 + source/getprm.f | 1 + source/gradient.f | 52 ++--- source/gradrot.f | 37 +-- source/hessian.f | 6 +- source/initatom.f | 1 + source/initprm.f | 4 +- source/initres.f | 1 + source/isokinetic.f | 476 ++++++++++++++++++++++++++++++++++++++ source/kangang.f | 7 +- source/kangle.f | 12 +- source/kangtor.f | 4 +- source/kbond.f | 4 +- source/kcharge.f | 10 +- source/kdipole.f | 6 +- source/keys.f | 5 +- source/kimprop.f | 7 +- source/kimptor.f | 9 +- source/kmpole.f | 35 ++- source/kopbend.f | 6 +- source/kopdist.f | 6 +- source/kpitors.f | 5 +- source/kpolar.f | 30 +-- source/ksolv.f | 50 +--- source/kstrbnd.f | 4 +- source/kstrtor.f | 4 +- source/ktors.f | 12 +- source/ktortor.f | 2 +- source/kurey.f | 6 +- source/kvdw.f | 8 +- source/lbfgs.f | 69 +++--- source/lights.f | 51 ++--- source/mdinit.f | 131 +++++++++-- source/mdstuf.f | 3 +- source/moldyn.f | 7 +- source/nblist.f | 542 ++++++++++---------------------------------- source/nucleic.f | 1 + source/params.f | 5 +- source/protein.f | 1 + source/ptable.f | 5 +- source/readgau.f | 5 +- source/readpdb.f | 3 + source/readseq.f | 1 + source/resdue.f | 94 ++++---- source/rgddyn.f | 19 +- source/rings.f | 16 +- source/scales.f | 2 +- source/shakeup.f | 1 + source/sizes.f | 34 ++- source/solute.f | 41 ++-- source/stochisok.f | 486 +++++++++++++++++++++++++++++++++++++++ source/testgrad.f | 2 +- source/testhess.f | 2 +- source/testpair.f | 12 +- source/testrot.f | 2 +- source/tncg.f | 88 +++---- source/torsions.f | 3 - source/vibbig.f | 6 +- source/warp.f | 5 +- 80 files changed, 1892 insertions(+), 1392 deletions(-) create mode 100644 source/isokinetic.f create mode 100644 source/stochisok.f diff --git a/source/analysis.f b/source/analysis.f index ca18f5e36..880386037 100644 --- a/source/analysis.f +++ b/source/analysis.f @@ -70,32 +70,32 @@ subroutine analysis (energy) c if (first) then first = .false. - if (.not. allocated(aesum)) allocate (aesum(n)) - if (.not. allocated(aeb)) allocate (aeb(n)) - if (.not. allocated(aea)) allocate (aea(n)) - if (.not. allocated(aeba)) allocate (aeba(n)) - if (.not. allocated(aeub)) allocate (aeub(n)) - if (.not. allocated(aeaa)) allocate (aeaa(n)) - if (.not. allocated(aeopb)) allocate (aeopb(n)) - if (.not. allocated(aeopd)) allocate (aeopd(n)) - if (.not. allocated(aeid)) allocate (aeid(n)) - if (.not. allocated(aeit)) allocate (aeit(n)) - if (.not. allocated(aet)) allocate (aet(n)) - if (.not. allocated(aept)) allocate (aept(n)) - if (.not. allocated(aebt)) allocate (aebt(n)) - if (.not. allocated(aeat)) allocate (aeat(n)) - if (.not. allocated(aett)) allocate (aett(n)) - if (.not. allocated(aev)) allocate (aev(n)) - if (.not. allocated(aec)) allocate (aec(n)) - if (.not. allocated(aecd)) allocate (aecd(n)) - if (.not. allocated(aed)) allocate (aed(n)) - if (.not. allocated(aem)) allocate (aem(n)) - if (.not. allocated(aep)) allocate (aep(n)) - if (.not. allocated(aer)) allocate (aer(n)) - if (.not. allocated(aes)) allocate (aes(n)) - if (.not. allocated(aelf)) allocate (aelf(n)) - if (.not. allocated(aeg)) allocate (aeg(n)) - if (.not. allocated(aex)) allocate (aex(n)) + if (.not. allocated(aesum)) allocate (aesum(maxatm)) + if (.not. allocated(aeb)) allocate (aeb(maxatm)) + if (.not. allocated(aea)) allocate (aea(maxatm)) + if (.not. allocated(aeba)) allocate (aeba(maxatm)) + if (.not. allocated(aeub)) allocate (aeub(maxatm)) + if (.not. allocated(aeaa)) allocate (aeaa(maxatm)) + if (.not. allocated(aeopb)) allocate (aeopb(maxatm)) + if (.not. allocated(aeopd)) allocate (aeopd(maxatm)) + if (.not. allocated(aeid)) allocate (aeid(maxatm)) + if (.not. allocated(aeit)) allocate (aeit(maxatm)) + if (.not. allocated(aet)) allocate (aet(maxatm)) + if (.not. allocated(aept)) allocate (aept(maxatm)) + if (.not. allocated(aebt)) allocate (aebt(maxatm)) + if (.not. allocated(aeat)) allocate (aeat(maxatm)) + if (.not. allocated(aett)) allocate (aett(maxatm)) + if (.not. allocated(aev)) allocate (aev(maxatm)) + if (.not. allocated(aec)) allocate (aec(maxatm)) + if (.not. allocated(aecd)) allocate (aecd(maxatm)) + if (.not. allocated(aed)) allocate (aed(maxatm)) + if (.not. allocated(aem)) allocate (aem(maxatm)) + if (.not. allocated(aep)) allocate (aep(maxatm)) + if (.not. allocated(aer)) allocate (aer(maxatm)) + if (.not. allocated(aes)) allocate (aes(maxatm)) + if (.not. allocated(aelf)) allocate (aelf(maxatm)) + if (.not. allocated(aeg)) allocate (aeg(maxatm)) + if (.not. allocated(aex)) allocate (aex(maxatm)) end if c c zero out energy partitioning components for each atom diff --git a/source/angang.f b/source/angang.f index c9e8418b1..75533ccc5 100644 --- a/source/angang.f +++ b/source/angang.f @@ -18,6 +18,7 @@ c c module angang + use sizes implicit none integer nangang integer, allocatable :: iaa(:,:) diff --git a/source/angbnd.f b/source/angbnd.f index f6f12bb07..37d787f1c 100644 --- a/source/angbnd.f +++ b/source/angbnd.f @@ -20,6 +20,7 @@ c c module angbnd + use sizes implicit none integer nangle integer, allocatable :: iang(:,:) diff --git a/source/angles.f b/source/angles.f index 1f38391ef..6c13a4a44 100644 --- a/source/angles.f +++ b/source/angles.f @@ -27,15 +27,13 @@ subroutine angles use iounit implicit none integer i,j,k,m - integer maxang c c c perform dynamic allocation of some global arrays c - maxang = 6 * n if (.not. allocated(iang)) allocate (iang(4,maxang)) if (.not. allocated(anglist)) - & allocate (anglist(maxval*(maxval-1)/2,n)) + & allocate (anglist(maxval*(maxval-1)/2,maxatm)) c c loop over all atoms, storing the atoms in each bond angle c diff --git a/source/attach.f b/source/attach.f index e0dd2faf7..21c10e965 100644 --- a/source/attach.f +++ b/source/attach.f @@ -35,12 +35,12 @@ subroutine attach maxn13 = 3 * maxval maxn14 = 9 * maxval maxn15 = 27 * maxval - if (.not. allocated(n13)) allocate (n13(n)) - if (.not. allocated(n14)) allocate (n14(n)) - if (.not. allocated(n15)) allocate (n15(n)) - if (.not. allocated(i13)) allocate (i13(maxn13,n)) - if (.not. allocated(i14)) allocate (i14(maxn14,n)) - if (.not. allocated(i15)) allocate (i15(maxn15,n)) + if (.not. allocated(n13)) allocate (n13(maxatm)) + if (.not. allocated(n14)) allocate (n14(maxatm)) + if (.not. allocated(n15)) allocate (n15(maxatm)) + if (.not. allocated(i13)) allocate (i13(maxn13,maxatm)) + if (.not. allocated(i14)) allocate (i14(maxn14,maxatm)) + if (.not. allocated(i15)) allocate (i15(maxn15,maxatm)) c c loop over all atoms finding all the 1-3 relationships; c note "n12" and "i12" have already been setup elsewhere diff --git a/source/bar.f b/source/bar.f index 8a1c9651e..2e3e57510 100644 --- a/source/bar.f +++ b/source/bar.f @@ -16,8 +16,8 @@ c via the Zwanzig free energy perturbation (FEP) and Bennett c acceptance ratio (BAR) methods c -c note current version takes as input a trajectory archive and -c key file for state A, and similar for state B; then finds the +c current version takes as input the trajectory archives and key +c files for state A, and similar for state B; then computes the c total potential energy for all frames of each trajectory under c control of both key files; finally the FEP and BAR algorithms c are used to compute the free energy for state A --> state B diff --git a/source/bath.f b/source/bath.f index a3e43216d..f7c8aba42 100644 --- a/source/bath.f +++ b/source/bath.f @@ -35,18 +35,31 @@ c thermostat choice of temperature control method to be used c barostat choice of pressure control method to be used c volscale choice of scaling method for Monte Carlo barostat +c q_iso1/2 masses of isokinetic thermostats +c v_iso1/2 velocities of isokinetic thermostats +c f_iso1/2 forces acting on isokinetic thermostats +c len_nhc In isokinetic dynamics, this is the Nose-Hoover chain length M. +c In stochastic isokinetic dynamics, this is the chain length L c c module bath implicit none integer maxnose parameter (maxnose=4) + integer len_nhc,isok_L,isok_M + integer isok_chain + parameter (isok_chain=4) integer voltrial real*8 kelvin,atmsph real*8 tautemp,taupres real*8 compress,collide real*8 eta,volmove real*8 vbar,qbar,gbar + real*8 q_iso1(isok_chain,3,10000) + real*8 q_iso2(isok_chain,3,10000) + real*8 v_iso1(isok_chain,3,10000) + real*8 v_iso2(isok_chain,3,10000) + real*8 f_iso1(isok_chain,3,10000) real*8 vnh(maxnose) real*8 qnh(maxnose) real*8 gnh(maxnose) diff --git a/source/bitors.f b/source/bitors.f index c5d02f1c4..56ee9260d 100644 --- a/source/bitors.f +++ b/source/bitors.f @@ -20,19 +20,16 @@ subroutine bitors use sizes use angbnd - use atoms use bitor use couple use iounit implicit none integer i,j,k integer ia,ib,ic,id,ie - integer maxbitor c c c perform dynamic allocation of some global arrays c - maxbitor = 8 * n if (.not. allocated(ibitor)) allocate (ibitor(5,maxbitor)) c c loop over all angles, storing the atoms in each bitorsion diff --git a/source/bonds.f b/source/bonds.f index ad91307c1..48d0f5e3e 100644 --- a/source/bonds.f +++ b/source/bonds.f @@ -25,14 +25,12 @@ subroutine bonds use iounit implicit none integer i,j,k,m - integer maxbnd c c c perform dynamic allocation of some global arrays c - maxbnd = 4 * n if (.not. allocated(ibnd)) allocate (ibnd(2,maxbnd)) - if (.not. allocated(bndlist)) allocate (bndlist(maxval,n)) + if (.not. allocated(bndlist)) allocate (bndlist(maxval,maxatm)) c c loop over all atoms, storing the atoms in each bond c diff --git a/source/cluster.f b/source/cluster.f index 202d7a21b..c0f967782 100644 --- a/source/cluster.f +++ b/source/cluster.f @@ -238,7 +238,7 @@ subroutine cluster c output the final list of atoms in each group c if (debug .and. use_group) then - do i = 1, ngrp + do i = 0, ngrp size = igrp(2,i) - igrp(1,i) + 1 if (size .ne. 0) then write (iout,50) i diff --git a/source/connolly.f b/source/connolly.f index afe6dafab..6d5600b48 100644 --- a/source/connolly.f +++ b/source/connolly.f @@ -44,7 +44,7 @@ c nfp number of convex faces c cynep number of convex edges in cycle c -c axyz atomic coordinates +c a atomic coordinates c ar atomic radii c pr probe radius c skip if true, atom is not used @@ -72,7 +72,7 @@ c c p probe coordinates c pa probe atom numbers -c vxyz vertex coordinates +c v vertex coordinates c va vertex atom number c vp vertex probe number c c circle center @@ -105,70 +105,6 @@ subroutine connolly (volume,area,radius,probe,exclude) real*8 radius(*) c c -c dimensions for arrays used by Connolly routines -c - maxcls = 50 * n - maxtt = 25 * n - maxt = 3 * n - maxp = 2 * n - maxv = 5 * n - maxen = 5 * n - maxfn = 2 * n - maxc = 5 * n - maxep = 5 * n - maxfs = 3 * n - maxcy = n - mxcyep = 30 - maxfp = n - mxfpcy = 10 -c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(ar)) allocate (ar(n)) - if (.not. allocated(axyz)) allocate (axyz(3,n)) - if (.not. allocated(skip)) allocate (skip(n)) - if (.not. allocated(nosurf)) allocate (nosurf(n)) - if (.not. allocated(afree)) allocate (afree(n)) - if (.not. allocated(abur)) allocate (abur(n)) - if (.not. allocated(cls)) allocate (cls(maxcls)) - if (.not. allocated(clst)) allocate (clst(maxcls)) - if (.not. allocated(acls)) allocate (acls(2,n)) - if (.not. allocated(ttfe)) allocate (ttfe(maxtt)) - if (.not. allocated(ttle)) allocate (ttle(maxtt)) - if (.not. allocated(enext)) allocate (enext(maxen)) - if (.not. allocated(tta)) allocate (tta(2,maxtt)) - if (.not. allocated(ttbur)) allocate (ttbur(maxtt)) - if (.not. allocated(ttfree)) allocate (ttfree(maxtt)) - if (.not. allocated(tfe)) allocate (tfe(maxt)) - if (.not. allocated(ta)) allocate (ta(2,maxt)) - if (.not. allocated(tr)) allocate (tr(maxt)) - if (.not. allocated(t)) allocate (t(3,maxt)) - if (.not. allocated(tax)) allocate (tax(3,maxt)) - if (.not. allocated(tfree)) allocate (tfree(maxt)) - if (.not. allocated(pa)) allocate (pa(3,maxp)) - if (.not. allocated(p)) allocate (p(3,maxp)) - if (.not. allocated(va)) allocate (va(maxv)) - if (.not. allocated(vp)) allocate (vp(maxv)) - if (.not. allocated(vxyz)) allocate (vxyz(3,maxv)) - if (.not. allocated(env)) allocate (env(2,maxen)) - if (.not. allocated(fnen)) allocate (fnen(3,maxfn)) - if (.not. allocated(ca)) allocate (ca(maxc)) - if (.not. allocated(ct)) allocate (ct(maxc)) - if (.not. allocated(cr)) allocate (cr(maxc)) - if (.not. allocated(c)) allocate (c(3,maxc)) - if (.not. allocated(epc)) allocate (epc(maxep)) - if (.not. allocated(epv)) allocate (epv(2,maxep)) - if (.not. allocated(afe)) allocate (afe(n)) - if (.not. allocated(ale)) allocate (ale(n)) - if (.not. allocated(epnext)) allocate (epnext(maxep)) - if (.not. allocated(fsen)) allocate (fsen(2,maxfs)) - if (.not. allocated(fsep)) allocate (fsep(2,maxfs)) - if (.not. allocated(cynep)) allocate (cynep(maxcy)) - if (.not. allocated(cyep)) allocate (cyep(mxcyep,maxcy)) - if (.not. allocated(fpa)) allocate (fpa(maxfp)) - if (.not. allocated(fpncy)) allocate (fpncy(maxfp)) - if (.not. allocated(fpcy)) allocate (fpcy(mxfpcy,maxfp)) -c c set the probe radius and the number of atoms c pr = probe @@ -178,9 +114,9 @@ subroutine connolly (volume,area,radius,probe,exclude) c radius ("exclude") is added to atomic radii c do i = 1, na - axyz(1,i) = x(i) - axyz(2,i) = y(i) - axyz(3,i) = z(i) + a(1,i) = x(i) + a(2,i) = y(i) + a(3,i) = z(i) ar(i) = radius(i) if (ar(i) .eq. 0.0d0) then skip(i) = .true. @@ -216,6 +152,7 @@ subroutine connolly (volume,area,radius,probe,exclude) c c subroutine wiggle + use sizes use faces implicit none integer i @@ -228,9 +165,9 @@ subroutine wiggle size = 0.000001d0 do i = 1, na call ranvec (vector) - axyz(1,i) = axyz(1,i) + size*vector(1) - axyz(2,i) = axyz(2,i) + size*vector(2) - axyz(3,i) = axyz(3,i) + size*vector(3) + a(1,i) = a(1,i) + size*vector(1) + a(2,i) = a(2,i) + size*vector(2) + a(3,i) = a(3,i) + size*vector(3) end do return end @@ -258,10 +195,12 @@ subroutine wiggle c c subroutine nearby + use sizes use faces implicit none - integer maxclsa + integer maxclsa,maxcube parameter (maxclsa=1000) + parameter (maxcube=40) integer i,j,k,m integer iptr,juse integer i1,j1,k1 @@ -271,19 +210,18 @@ subroutine nearby integer jcls,jmin integer jmincls,jmold integer ncls,nclsa - integer maxcube integer clsa(maxclsa) integer itnl(maxclsa) integer, allocatable :: icuptr(:) integer, allocatable :: ico(:,:) - integer, allocatable :: icube(:,:,:) + integer icube(maxcube,maxcube,maxcube) real*8 radmax,width real*8 sum,sumi real*8 dist2,d2,r2 real*8 vect1,vect2,vect3 real*8 comin(3) - logical, allocatable :: scube(:,:,:) - logical, allocatable :: sscube(:,:,:) + logical scube(maxcube,maxcube,maxcube) + logical sscube(maxcube,maxcube,maxcube) c c c ignore all atoms that are completely inside another atom; @@ -292,7 +230,7 @@ subroutine nearby do i = 1, na-1 if (.not. skip(i)) then do j = i+1, na - d2 = dist2(axyz(1,i),axyz(1,j)) + d2 = dist2(a(1,i),a(1,j)) r2 = (ar(i) - ar(j))**2 if (.not.skip(j) .and. d2.lt.r2) then if (ar(i) .lt. ar(j)) then @@ -309,11 +247,11 @@ subroutine nearby c radmax = 0.0d0 do k = 1, 3 - comin(k) = axyz(k,1) + comin(k) = a(k,1) end do do i = 1, na do k = 1, 3 - if (axyz(k,i) .lt. comin(k)) comin(k) = axyz(k,i) + if (a(k,i) .lt. comin(k)) comin(k) = a(k,i) end do if (ar(i) .gt. radmax) radmax = ar(i) end do @@ -325,18 +263,14 @@ subroutine nearby c c perform dynamic allocation of some local arrays c - maxcube = 40 allocate (icuptr(na)) allocate (ico(3,na)) - allocate (icube(maxcube,maxcube,maxcube)) - allocate (scube(maxcube,maxcube,maxcube)) - allocate (sscube(maxcube,maxcube,maxcube)) c c set up cube arrays; first the integer coordinate arrays c do i = 1, na do k = 1, 3 - ico(k,i) = (axyz(k,i)-comin(k))/width + 1 + ico(k,i) = (a(k,i)-comin(k))/width + 1 if (ico(k,i) .lt. 1) then call cerror ('Cube Coordinate Too Small') else if (ico(k,i) .gt. maxcube) then @@ -387,7 +321,7 @@ subroutine nearby c c check for duplicate atoms, turn off one of them c - if (dist2(axyz(1,iatom),axyz(1,iptr)) .le. 0.0d0) then + if (dist2(a(1,iatom),a(1,iptr)) .le. 0.0d0) then skip(iatom) = .true. goto 30 end if @@ -466,11 +400,11 @@ subroutine nearby c distance check c sum = sumi + ar(j) - vect1 = abs(axyz(1,j) - axyz(1,i)) + vect1 = abs(a(1,j) - a(1,i)) if (vect1 .ge. sum) goto 50 - vect2 = abs(axyz(2,j) - axyz(2,i)) + vect2 = abs(a(2,j) - a(2,i)) if (vect2 .ge. sum) goto 50 - vect3 = abs(axyz(3,j) - axyz(3,i)) + vect3 = abs(a(3,j) - a(3,i)) if (vect3 .ge. sum) goto 50 d2 = vect1**2 + vect2**2 + vect3**2 if (d2 .ge. sum**2) goto 50 @@ -537,9 +471,6 @@ subroutine nearby c deallocate (icuptr) deallocate (ico) - deallocate (icube) - deallocate (scube) - deallocate (sscube) return end c @@ -556,6 +487,7 @@ subroutine nearby c c subroutine torus + use sizes use faces implicit none integer ia,ja,jn @@ -761,7 +693,7 @@ subroutine place call gettor (ia,ja,ttok,bij,hij,uij) do km = 1, nmnb ka = mnb(km) - discls(km) = dist2(bij,axyz(1,ka)) + discls(km) = dist2(bij,a(1,ka)) sumcls(km) = (pr+ar(ka))**2 c c initialize link to next farthest out neighbor @@ -870,7 +802,7 @@ subroutine place c c compare distance to sum of radii c - d2 = dist2(pijk,axyz(1,la)) + d2 = dist2(pijk,a(1,la)) if (d2 .le. sumcls(lm)) goto 110 90 continue lm = lkcls(lm) @@ -900,26 +832,26 @@ subroutine place c if (nv+3 .gt. maxv) call cerror ('Too many Vertices') do k = 1, 3 - vxyz(k,nv+1) = axyz(k,ia) - p(k,np) - vxyz(k,nv+2) = axyz(k,ja) - p(k,np) - vxyz(k,nv+3) = axyz(k,ka) - p(k,np) + v(k,nv+1) = a(k,ia) - p(k,np) + v(k,nv+2) = a(k,ja) - p(k,np) + v(k,nv+3) = a(k,ka) - p(k,np) end do c c calculate determinant of vectors defining triangle c - det = vxyz(1,nv+1)*vxyz(2,nv+2)*vxyz(3,nv+3) - & + vxyz(1,nv+2)*vxyz(2,nv+3)*vxyz(3,nv+1) - & + vxyz(1,nv+3)*vxyz(2,nv+1)*vxyz(3,nv+2) - & - vxyz(1,nv+3)*vxyz(2,nv+2)*vxyz(3,nv+1) - & - vxyz(1,nv+2)*vxyz(2,nv+1)*vxyz(3,nv+3) - & - vxyz(1,nv+1)*vxyz(2,nv+3)*vxyz(3,nv+2) + det = v(1,nv+1)*v(2,nv+2)*v(3,nv+3) + & + v(1,nv+2)*v(2,nv+3)*v(3,nv+1) + & + v(1,nv+3)*v(2,nv+1)*v(3,nv+2) + & - v(1,nv+3)*v(2,nv+2)*v(3,nv+1) + & - v(1,nv+2)*v(2,nv+1)*v(3,nv+3) + & - v(1,nv+1)*v(2,nv+3)*v(3,nv+2) c c now add probe coordinates to vertices c do k = 1, 3 - vxyz(k,nv+1) = p(k,np) + vxyz(k,nv+1)*pr/(ar(ia)+pr) - vxyz(k,nv+2) = p(k,np) + vxyz(k,nv+2)*pr/(ar(ja)+pr) - vxyz(k,nv+3) = p(k,np) + vxyz(k,nv+3)*pr/(ar(ka)+pr) + v(k,nv+1) = p(k,np) + v(k,nv+1)*pr/(ar(ia)+pr) + v(k,nv+2) = p(k,np) + v(k,nv+2)*pr/(ar(ja)+pr) + v(k,nv+3) = p(k,np) + v(k,nv+3)*pr/(ar(ka)+pr) end do c c want the concave face to have counter-clockwise orientation @@ -929,9 +861,9 @@ subroutine place c swap second and third vertices c do k = 1, 3 - tempv(k) = vxyz(k,nv+2) - vxyz(k,nv+2) = vxyz(k,nv+3) - vxyz(k,nv+3) = tempv(k) + tempv(k) = v(k,nv+2) + v(k,nv+2) = v(k,nv+3) + v(k,nv+3) = tempv(k) end do c c set up pointers from probe to atoms @@ -1184,7 +1116,7 @@ subroutine saddles c vector from atom to torus center c do k = 1, 3 - atvect(k) = t(k,it) - axyz(k,ia) + atvect(k) = t(k,it) - a(k,ia) end do factor = ar(ia) / (ar(ia)+pr) c @@ -1196,12 +1128,15 @@ subroutine saddles c circle center c do k = 1, 3 - c(k,nc) = axyz(k,ia) + factor*atvect(k) + c(k,nc) = a(k,ia) + factor*atvect(k) end do c -c pointer from circle to atom and to torus +c pointer from circle to atom c ca(nc) = ia +c +c pointer from circle to torus +c ct(nc) = it c c circle radius @@ -1509,7 +1444,7 @@ subroutine saddles c fsep(2,nfs) = nep c -c nothing to do for buried torus +c buried torus; do nothing with it c 80 continue end do @@ -1533,9 +1468,8 @@ subroutine gettor (ia,ja,ttok,torcen,torad,torax) use faces implicit none integer k,ia,ja - real*8 dist2,dij - real*8 temp,temp1 - real*8 temp2 + real*8 dist2,dij,temp + real*8 temp1,temp2 real*8 torad real*8 torcen(3) real*8 torax(3) @@ -1548,12 +1482,12 @@ subroutine gettor (ia,ja,ttok,torcen,torad,torax) c get the distance between the two atoms c ttok = .false. - dij = sqrt(dist2(axyz(1,ia),axyz(1,ja))) + dij = sqrt(dist2(a(1,ia),a(1,ja))) c c find a unit vector along interatomic (torus) axis c do k = 1, 3 - vij(k) = axyz(k,ja) - axyz(k,ia) + vij(k) = a(k,ja) - a(k,ia) uij(k) = vij(k) / dij end do c @@ -1561,7 +1495,7 @@ subroutine gettor (ia,ja,ttok,torcen,torad,torax) c temp = 1.0d0 + ((ar(ia)+pr)**2-(ar(ja)+pr)**2)/dij**2 do k = 1, 3 - bij(k) = axyz(k,ia) + 0.5d0*vij(k)*temp + bij(k) = a(k,ia) + 0.5d0*vij(k)*temp end do c c skip if atoms too far apart (should not happen) @@ -1623,7 +1557,7 @@ subroutine getprb (ia,ja,ka,prbok,tb,bijk,hijk,uijk) tb = .false. call gettor (ia,ja,tok,tij,rij,uij) if (.not. tok) return - dat2 = dist2(axyz(1,ka),tij) + dat2 = dist2(a(1,ka),tij) rad2 = (ar(ka)+pr)**2 - rij**2 c c if "ka" less than "ja", then all we care about @@ -1662,7 +1596,7 @@ subroutine getprb (ia,ja,ka,prbok,tb,bijk,hijk,uijk) do k = 1, 3 bijk(k) = tij(k) + utb(k)*fact end do - dba = dist2(axyz(1,ia),bijk) + dba = dist2(a(1,ia),bijk) rip2 = (ar(ia) + pr)**2 rad = rip2 - dba if (rad .lt. 0.0d0) then @@ -1736,10 +1670,8 @@ subroutine contact integer ncypa,icya,jcya,kcya integer ncyep,icyep,jcyep integer ncyold,nused,lookv - integer aic(maxepa) - integer aia(maxepa) - integer aep(maxepa) - integer av(2,maxepa) + integer aic(maxepa),aia(maxepa) + integer aep(maxepa),av(2,maxepa) integer ncyepa(maxcypa) integer cyepa(mxcyep,maxcypa) real*8 anorm,anaa,factor @@ -1830,8 +1762,8 @@ subroutine contact c sometimes we use one vector, sometimes the other c do k = 1, 3 - acvect(k,nepa) = c(k,ic) - axyz(k,ia) - aavect(k,nepa) = axyz(k,ia2) - axyz(k,ia) + acvect(k,nepa) = c(k,ic) - a(k,ia) + aavect(k,nepa) = a(k,ia2) - a(k,ia) end do c c circle radius @@ -1966,7 +1898,8 @@ subroutine contact end if 70 continue c -c this cycle is finished, store number of edges in cycle +c this cycle is finished +c store number of edges in cycle c ncyepa(ncypa) = ncyep cynep(ncy) = ncyep @@ -1983,7 +1916,13 @@ subroutine contact do icya = 1, ncypa do jcya = 1, ncypa jcy = ncyold + jcya +c +c initialize +c cycy(icya,jcya) = .true. +c +c check for same cycle +c if (icya .eq. jcya) goto 90 c c if cycle j has two or fewer edges, nothing can @@ -2013,7 +1952,7 @@ subroutine contact c north pole and unit vector pointing south c do k = 1, 3 - pole(k) = factor*aavect(k,iepa) + axyz(k,ia) + pole(k) = factor*aavect(k,iepa) + a(k,ia) unvect(k) = -aavect(k,iepa) / anaa end do cycy(icya,jcya) = ptincy(pole,unvect,jcy) @@ -2366,7 +2305,7 @@ subroutine vam (volume,area) end do do ke = 1, 3 do k = 1, 3 - vects(k,ke) = vxyz(k,ivs(ke)) - p(k,ip) + vects(k,ke) = v(k,ivs(ke)) - p(k,ip) end do end do c @@ -2500,8 +2439,8 @@ subroutine vam (volume,area) iv1 = env(1,ien) iv2 = env(2,ien) do k = 1, 3 - vect3(k) = vxyz(k,iv1) - fncen(k,ifn) - vect4(k) = vxyz(k,iv2) - fncen(k,ifn) + vect3(k) = v(k,iv1) - fncen(k,ifn) + vect4(k) = v(k,iv2) - fncen(k,ifn) end do do ke2 = 1, 3 if (ispind(ke) .eq. ispnd2(ke2)) goto 40 @@ -2514,8 +2453,8 @@ subroutine vam (volume,area) iv1 = env(1,ien) iv2 = env(2,ien) do k = 1, 3 - vect7(k) = vxyz(k,iv1) - fncen(k,jfn) - vect8(k) = vxyz(k,iv2) - fncen(k,jfn) + vect7(k) = v(k,iv1) - fncen(k,jfn) + vect8(k) = v(k,iv2) - fncen(k,jfn) end do c c check whether point lies on spindle arc @@ -2554,8 +2493,8 @@ subroutine vam (volume,area) iv1 = env(1,ien) iv2 = env(2,ien) do k = 1, 3 - vect3(k) = vxyz(k,iv1) - fncen(k,ifn) - vect4(k) = vxyz(k,iv2) - fncen(k,ifn) + vect3(k) = v(k,iv1) - fncen(k,ifn) + vect4(k) = v(k,iv2) - fncen(k,ifn) end do do ke2 = 1, 3 if (ispind(ke) .eq. ispnd2(ke2)) goto 70 @@ -2568,8 +2507,8 @@ subroutine vam (volume,area) iv1 = env(1,ien) iv2 = env(2,ien) do k = 1, 3 - vect7(k) = vxyz(k,iv1) - fncen(k,jfn) - vect8(k) = vxyz(k,iv2) - fncen(k,jfn) + vect7(k) = v(k,iv1) - fncen(k,jfn) + vect8(k) = v(k,iv2) - fncen(k,jfn) end do c c check whether point lies on spindle arc @@ -2630,11 +2569,11 @@ subroutine vam (volume,area) ien = fnen(kv,jfn) iv = env(1,ien) do k = 1, 3 - vect1(k) = vxyz(k,iv) - fncen(k,ifn) + vect1(k) = v(k,iv) - fncen(k,ifn) end do call vnorm (vect1,vect1) do ke = 1, 3 - dt = dot(fnvect(1,ke,ifn),vxyz(1,iv)) + dt = dot(fnvect(1,ke,ifn),v(1,iv)) if (dt .gt. 0.0d0) goto 80 end do vip(kv) = .true. @@ -2831,6 +2770,7 @@ subroutine vam (volume,area) c c function depth (ip,alt) + use sizes use faces implicit none integer k,ip,ia1,ia2,ia3 @@ -2843,9 +2783,9 @@ function depth (ip,alt) ia2 = pa(2,ip) ia3 = pa(3,ip) do k = 1, 3 - vect1(k) = axyz(k,ia1) - axyz(k,ia3) - vect2(k) = axyz(k,ia2) - axyz(k,ia3) - vect3(k) = p(k,ip) - axyz(k,ia3) + vect1(k) = a(k,ia1) - a(k,ia3) + vect2(k) = a(k,ia2) - a(k,ia3) + vect3(k) = p(k,ip) - a(k,ia3) end do call vcross (vect1,vect2,vect4) call vnorm (vect4,vect4) @@ -2869,15 +2809,12 @@ function depth (ip,alt) c c subroutine measpm (ifn,prism) + use sizes use faces implicit none - integer k,ke,ien - integer iv,ia,ip,ifn - real*8 prism,height - real*8 vect1(3) - real*8 vect2(3) - real*8 vect3(3) - real*8 pav(3,3) + integer k,ke,ien,iv,ia,ip,ifn + real*8 prism,height,pav(3,3) + real*8 vect1(3),vect2(3),vect3(3) c c height = 0.0d0 @@ -2885,10 +2822,10 @@ subroutine measpm (ifn,prism) ien = fnen(ke,ifn) iv = env(1,ien) ia = va(iv) - height = height + axyz(3,ia) + height = height + a(3,ia) ip = vp(iv) do k = 1, 3 - pav(k,ke) = axyz(k,ia) - p(k,ip) + pav(k,ke) = a(k,ia) - p(k,ip) end do end do height = height / 3.0d0 @@ -2910,6 +2847,7 @@ subroutine measpm (ifn,prism) c c subroutine measfp (ifp,areap,volp) + use sizes use faces use math implicit none @@ -2917,18 +2855,15 @@ subroutine measfp (ifp,areap,volp) integer ia,ia2,ic integer it,iv1,iv2 integer ncycle,ieuler - integer icyptr,icy - integer nedge + integer icyptr,icy,nedge real*8 areap,volp real*8 dot,dt,gauss real*8 vecang,angle,geo real*8 pcurve,gcurve - real*8 vect1(3) - real*8 vect2(3) - real*8 acvect(3) - real*8 aavect(3) - real*8 radial(3,mxcyep) + real*8 vect1(3),vect2(3) + real*8 acvect(3),aavect(3) real*8 tanv(3,2,mxcyep) + real*8 radial(3,mxcyep) c c ia = fpa(ifp) @@ -2953,8 +2888,8 @@ subroutine measfp (ifp,areap,volp) ia2 = ta(1,it) end if do k = 1, 3 - acvect(k) = c(k,ic) - axyz(k,ia) - aavect(k) = axyz(k,ia2) - axyz(k,ia) + acvect(k) = c(k,ic) - a(k,ia) + aavect(k) = a(k,ia2) - a(k,ia) end do call vnorm (aavect,aavect) dt = dot(acvect,aavect) @@ -2965,9 +2900,9 @@ subroutine measfp (ifp,areap,volp) angle = 2.0d0 * pi else do k = 1, 3 - vect1(k) = vxyz(k,iv1) - c(k,ic) - vect2(k) = vxyz(k,iv2) - c(k,ic) - radial(k,ke) = vxyz(k,iv1) - axyz(k,ia) + vect1(k) = v(k,iv1) - c(k,ic) + vect2(k) = v(k,iv2) - c(k,ic) + radial(k,ke) = v(k,iv1) - a(k,ia) end do call vnorm (radial(1,ke),radial(1,ke)) call vcross (vect1,aavect,tanv(1,1,ke)) @@ -3012,6 +2947,7 @@ subroutine measfp (ifp,areap,volp) c c subroutine measfs (ifs,areas,vols,areasp,volsp) + use sizes use faces use math implicit none @@ -3026,8 +2962,7 @@ subroutine measfs (ifs,areas,vols,areasp,volsp) real*8 theta1,theta2 real*8 rat,thetaq real*8 cone1,cone2 - real*8 term1,term2 - real*8 term3 + real*8 term1,term2,term3 real*8 spin,volt real*8 vect1(3) real*8 vect2(3) @@ -3041,7 +2976,7 @@ subroutine measfs (ifs,areas,vols,areasp,volsp) ia1 = ta(1,it) ia2 = ta(2,it) do k = 1, 3 - aavect(k) = axyz(k,ia2) - axyz(k,ia1) + aavect(k) = a(k,ia2) - a(k,ia1) end do call vnorm (aavect,aavect) iv1 = epv(1,iep) @@ -3050,14 +2985,14 @@ subroutine measfs (ifs,areas,vols,areasp,volsp) phi = 2.0d0 * pi else do k = 1, 3 - vect1(k) = vxyz(k,iv1) - c(k,ic) - vect2(k) = vxyz(k,iv2) - c(k,ic) + vect1(k) = v(k,iv1) - c(k,ic) + vect2(k) = v(k,iv2) - c(k,ic) end do phi = vecang(vect1,vect2,aavect,1.0d0) end if do k = 1, 3 - vect1(k) = axyz(k,ia1) - t(k,it) - vect2(k) = axyz(k,ia2) - t(k,it) + vect1(k) = a(k,ia1) - t(k,it) + vect2(k) = a(k,ia2) - t(k,it) end do d1 = -dot(vect1,aavect) d2 = dot(vect2,aavect) @@ -3095,8 +3030,8 @@ subroutine measfs (ifs,areas,vols,areasp,volsp) call cerror ('IA1 Inconsistency in MEASFS') end if do k = 1, 3 - vect1(k) = c(k,ic1) - axyz(k,ia1) - vect2(k) = c(k,ic2) - axyz(k,ia2) + vect1(k) = c(k,ic1) - a(k,ia1) + vect2(k) = c(k,ic2) - a(k,ia2) end do w1 = dot(vect1,aavect) w2 = -dot(vect2,aavect) @@ -3104,10 +3039,10 @@ subroutine measfs (ifs,areas,vols,areasp,volsp) cone2 = phi * (w2*cr(ic2)**2)/6.0d0 term1 = (tr(it)**2) * pr * (sin(theta1)+sin(theta2)) term2 = sin(theta1)*cos(theta1) + theta1 - & + sin(theta2)*cos(theta2) + theta2 + & + sin(theta2)*cos(theta2) + theta2 term2 = tr(it) * (pr**2) * term2 term3 = sin(theta1)*cos(theta1)**2 + 2.0d0*sin(theta1) - & + sin(theta2)*cos(theta2)**2 + 2.0d0*sin(theta2) + & + sin(theta2)*cos(theta2)**2 + 2.0d0*sin(theta2) term3 = (pr**3 / 3.0d0) * term3 volt = (phi/2.0d0) * (term1-term2+term3) vols = volt + cone1 + cone2 @@ -3131,6 +3066,7 @@ subroutine measfs (ifs,areas,vols,areasp,volsp) c c subroutine measfn (ifn,arean,voln) + use sizes use faces use math implicit none @@ -3152,8 +3088,8 @@ subroutine measfn (ifn,arean,voln) ia = va(iv) ip = vp(iv) do k = 1, 3 - pvv(k,ke) = vxyz(k,iv) - p(k,ip) - pav(k,ke) = axyz(k,ia) - p(k,ip) + pvv(k,ke) = v(k,iv) - p(k,ip) + pav(k,ke) = a(k,ia) - p(k,ip) end do if (pr .gt. 0.0d0) call vnorm (pvv(1,ke),pvv(1,ke)) end do @@ -3192,6 +3128,7 @@ subroutine measfn (ifn,arean,voln) c c subroutine projct (pnt,unvect,icy,ia,spv,nedge,fail) + use sizes use faces implicit none integer k,ke,icy,ia @@ -3217,7 +3154,7 @@ subroutine projct (pnt,unvect,icy,ia,spv,nedge,fail) c vector from north pole to vertex c do k = 1, 3 - polev(k) = vxyz(k,iv) - pnt(k) + polev(k) = v(k,iv) - pnt(k) end do c c calculate multiplication factor @@ -3253,13 +3190,13 @@ subroutine projct (pnt,unvect,icy,ia,spv,nedge,fail) c c function ptincy (pnt,unvect,icy) + use sizes use faces implicit none integer k,ke,icy,iep integer ic,it,iatom integer iaoth,nedge - real*8 dot,rotang - real*8 totang + real*8 dot,rotang,totang real*8 unvect(3) real*8 pnt(3) real*8 acvect(3) @@ -3282,7 +3219,7 @@ function ptincy (pnt,unvect,icy) iaoth = ta(1,it) end if do k = 1, 3 - acvect(k) = axyz(k,iaoth) - axyz(k,iatom) + acvect(k) = a(k,iaoth) - a(k,iatom) cpvect(k) = pnt(k) - c(k,ic) end do if (dot(acvect,cpvect) .ge. 0.0d0) then @@ -3572,12 +3509,9 @@ function vecang (v1,v2,axis,hand) use math implicit none real*8 vecang,hand - real*8 angle,dt - real*8 a1,a2,a12 - real*8 anorm,dot - real*8 triple - real*8 v1(3),v2(3) - real*8 axis(3) + real*8 angle,a1,a2,a12,dt + real*8 anorm,dot,triple + real*8 v1(3),v2(3),axis(3) c c a1 = anorm(v1) @@ -3612,24 +3546,15 @@ subroutine cirpln (circen,cirrad,cirvec,plncen,plnvec, & cinsp,cintp,xpnt1,xpnt2) implicit none integer k - real*8 anorm,dot - real*8 dcp,dir - real*8 ratio,rlen - real*8 cirrad - real*8 circen(3) - real*8 cirvec(3) - real*8 plncen(3) - real*8 plnvec(3) - real*8 xpnt1(3) - real*8 xpnt2(3) - real*8 cpvect(3) - real*8 pnt1(3) - real*8 vect1(3) - real*8 vect2(3) - real*8 uvect1(3) - real*8 uvect2(3) - logical cinsp - logical cintp + real*8 anorm,dot,dcp,dir + real*8 ratio,rlen,cirrad + real*8 circen(3),cirvec(3) + real*8 plncen(3),plnvec(3) + real*8 xpnt1(3),xpnt2(3) + real*8 cpvect(3),pnt1(3) + real*8 vect1(3),vect2(3) + real*8 uvect1(3),uvect2(3) + logical cinsp,cintp c c do k = 1, 3 @@ -3682,16 +3607,10 @@ subroutine cirpln (circen,cirrad,cirvec,plncen,plnvec, subroutine gendot (ndots,dots,radius,xcenter,ycenter,zcenter) use math implicit none - integer i,j,k - integer ndots - integer nequat - integer nvert - integer nhoriz - real*8 fi,fj - real*8 x,y,z,xy - real*8 xcenter - real*8 ycenter - real*8 zcenter + integer i,j,k,ndots + integer nequat,nvert,nhoriz + real*8 fi,fj,x,y,z,xy + real*8 xcenter,ycenter,zcenter real*8 radius real*8 dots(3,*) c diff --git a/source/diffuse.f b/source/diffuse.f index e4d5d8053..1f8ca6417 100644 --- a/source/diffuse.f +++ b/source/diffuse.f @@ -147,23 +147,12 @@ program diffuse do i = 1, 20 list(i) = 0 end do - i = 0 - dowhile (exist) - call nextarg (string,exist) - if (exist) then - read (string,*,err=100,end=100) list(i+1) - i = i + 1 - end if - end do - 100 continue - if (i .eq. 0) then - write (iout,110) - 110 format (/,' Numbers of any Atoms to be Removed : ',$) - read (input,120) record - 120 format (a120) - read (record,*,err=130,end=130) (list(i),i=1,20) - 130 continue - end if + write (iout,100) + 100 format (/,' Numbers of any Atoms to be Removed : ',$) + read (input,110) record + 110 format (a120) + read (record,*,err=120,end=120) (list(i),i=1,20) + 120 continue i = 1 do while (list(i) .ne. 0) list(i) = max(-n,min(n,list(i))) @@ -205,8 +194,8 @@ program diffuse end if end do nmol = k - write (iout,140) nmol - 140 format (/,' Total Number of Molecules :',i16) + write (iout,130) nmol + 130 format (/,' Total Number of Molecules :',i16) c c count the number of coordinate frames in the archive file c @@ -221,8 +210,8 @@ program diffuse rewind (unit=iarc) stop = min(nframe,stop) nframe = (stop-start)/step + 1 - write (iout,150) nframe - 150 format (/,' Number of Coordinate Frames :',i14) + write (iout,140) nframe + 140 format (/,' Number of Coordinate Frames :',i14) c c perform dynamic allocation of some local arrays c @@ -236,8 +225,8 @@ program diffuse c c get the archived coordinates for each frame in turn c - write (iout,160) - 160 format (/,' Reading the Coordinates Archive File :',/) + write (iout,150) + 150 format (/,' Reading the Coordinates Archive File :',/) nframe = 0 iframe = start skip = start @@ -248,11 +237,11 @@ program diffuse iframe = iframe + step skip = step call readxyz (iarc) - if (n .eq. 0) goto 180 + if (n .eq. 0) goto 170 nframe = nframe + 1 if (mod(nframe,100) .eq. 0) then - write (iout,170) nframe - 170 format (4x,'Processing Coordinate Frame',i13) + write (iout,160) nframe + 160 format (4x,'Processing Coordinate Frame',i13) end if c c unfold each molecule to get its corrected center of mass @@ -290,11 +279,11 @@ program diffuse zcm(i,nframe) = zold + zr end do end do - 180 continue + 170 continue close (unit=iarc) if (mod(nframe,100) .ne. 0) then - write (iout,190) nframe - 190 format (4x,'Processing Coordinate Frame',i13) + write (iout,180) nframe + 180 format (4x,'Processing Coordinate Frame',i13) end if c c increment the squared displacements for each frame pair @@ -339,8 +328,8 @@ program diffuse c c estimate the diffusion constant via the Einstein relation c - write (iout,200) - 200 format (/,' Mean Squared Displacements and Self-Diffusion', + write (iout,190) + 190 format (/,' Mean Squared Displacements and Self-Diffusion', & ' Constant :', & //,5x,'Time Gap',6x,'X MSD',7x,'Y MSD',7x,'Z MSD', & 7x,'R MSD',4x,'Diff Const', @@ -353,8 +342,8 @@ program diffuse zvalue = zmsd(i) / 2.0d0 rvalue = (xmsd(i) + ymsd(i) + zmsd(i)) / 6.0d0 dvalue = rvalue / delta - write (iout,210) delta,xvalue,yvalue,zvalue,rvalue,dvalue - 210 format (f12.2,4f12.2,f12.4) + write (iout,200) delta,xvalue,yvalue,zvalue,rvalue,dvalue + 200 format (f12.2,4f12.2,f12.4) end do c c perform deallocation of some local arrays diff --git a/source/domega.f b/source/domega.f index 0a02a2e52..44f692390 100644 --- a/source/domega.f +++ b/source/domega.f @@ -41,32 +41,33 @@ c c module domega + use sizes implicit none - real*8, allocatable :: tesum(:) - real*8, allocatable :: teb(:) - real*8, allocatable :: tea(:) - real*8, allocatable :: teba(:) - real*8, allocatable :: teub(:) - real*8, allocatable :: teaa(:) - real*8, allocatable :: teopb(:) - real*8, allocatable :: teopd(:) - real*8, allocatable :: teid(:) - real*8, allocatable :: teit(:) - real*8, allocatable :: tet(:) - real*8, allocatable :: tept(:) - real*8, allocatable :: tebt(:) - real*8, allocatable :: teat(:) - real*8, allocatable :: tett(:) - real*8, allocatable :: tev(:) - real*8, allocatable :: tec(:) - real*8, allocatable :: tecd(:) - real*8, allocatable :: ted(:) - real*8, allocatable :: tem(:) - real*8, allocatable :: tep(:) - real*8, allocatable :: ter(:) - real*8, allocatable :: tes(:) - real*8, allocatable :: telf(:) - real*8, allocatable :: teg(:) - real*8, allocatable :: tex(:) + real*8 tesum(maxrot) + real*8 teb(maxrot) + real*8 tea(maxrot) + real*8 teba(maxrot) + real*8 teub(maxrot) + real*8 teaa(maxrot) + real*8 teopb(maxrot) + real*8 teopd(maxrot) + real*8 teid(maxrot) + real*8 teit(maxrot) + real*8 tet(maxrot) + real*8 tept(maxrot) + real*8 tebt(maxrot) + real*8 teat(maxrot) + real*8 tett(maxrot) + real*8 tev(maxrot) + real*8 tec(maxrot) + real*8 tecd(maxrot) + real*8 ted(maxrot) + real*8 tem(maxrot) + real*8 tep(maxrot) + real*8 ter(maxrot) + real*8 tes(maxrot) + real*8 telf(maxrot) + real*8 teg(maxrot) + real*8 tex(maxrot) save end diff --git a/source/dynamic.f b/source/dynamic.f index 13884cfbc..f7cf9b4c4 100644 --- a/source/dynamic.f +++ b/source/dynamic.f @@ -261,9 +261,17 @@ program dynamic write (iout,390) 390 format (/,' Molecular Dynamics Trajectory via', & ' r-RESPA MTS Algorithm') - else + else if (integrate .eq. 'STOCH-ISOK') then write (iout,400) 400 format (/,' Molecular Dynamics Trajectory via', + & ' Stochastic Isokinetic Algorithm') + else if (integrate .eq. 'ISOKINETIC') then + write (iout,410) + 410 format (/,' Molecular Dynamics Trajectory via', + & ' Isokinetic NHC Algorithm') + else + write (iout,420) + 420 format (/,' Molecular Dynamics Trajectory via', & ' Modified Beeman Algorithm') end if c @@ -284,6 +292,10 @@ program dynamic call rgdstep (istep,dt) else if (integrate .eq. 'RESPA') then call respa (istep,dt) + else if (integrate .eq. 'ISOKINETIC') then + call isokinetic (istep,dt) + else if (integrate .eq. 'STOCH-ISOK') then + call stochisok (istep,dt) else call beeman (istep,dt) end if diff --git a/source/eangang2.f b/source/eangang2.f index ff482d4a9..55462bdad 100644 --- a/source/eangang2.f +++ b/source/eangang2.f @@ -250,10 +250,6 @@ subroutine eangang2a (i) yie = y(ie) zie = z(ie) c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(deaa)) allocate (deaa(3,n)) -c c zero out the first derivative components c deaa(1,ia) = 0.0d0 diff --git a/source/embed.f b/source/embed.f index 90258085d..108881308 100644 --- a/source/embed.f +++ b/source/embed.f @@ -2955,9 +2955,9 @@ function vdwerr (derivs) real*8 dx,dy,dz,gx,gy,gz real*8 dstsq,blosq real*8 radi,radsq - real*8, allocatable :: xsort(:) - real*8, allocatable :: ysort(:) - real*8, allocatable :: zsort(:) + real*8 xsort(maxlight) + real*8 ysort(maxlight) + real*8 zsort(maxlight) real*8 derivs(3,*) c c @@ -2969,9 +2969,6 @@ function vdwerr (derivs) c perform dynamic allocation of some local arrays c allocate (skip(n)) - allocate (xsort(n)) - allocate (ysort(n)) - allocate (zsort(n)) c c transfer coordinates and zero out atoms to be skipped c @@ -3043,9 +3040,6 @@ function vdwerr (derivs) c perform deallocation of some local arrays c deallocate (skip) - deallocate (xsort) - deallocate (ysort) - deallocate (zsort) return end c diff --git a/source/eopbend2.f b/source/eopbend2.f index 07d1c2593..2095c13d1 100644 --- a/source/eopbend2.f +++ b/source/eopbend2.f @@ -230,10 +230,6 @@ subroutine eopbend2a (i) yid = y(id) zid = z(id) c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(deopb)) allocate (deopb(3,n)) -c c zero out the first derivative components c deopb(1,ia) = 0.0d0 diff --git a/source/esolv1.f b/source/esolv1.f index 80fe89cdb..f3afcc78b 100644 --- a/source/esolv1.f +++ b/source/esolv1.f @@ -47,15 +47,11 @@ subroutine esolv1 eace = 0.0d0 do i = 1, n drb(i) = 0.0d0 + drbp(i) = 0.0d0 des(1,i) = 0.0d0 des(2,i) = 0.0d0 des(3,i) = 0.0d0 end do - if (solvtyp(1:2) .eq. 'GK') then - do i = 1, n - drbp(i) = 0.0d0 - end do - end if c c set a value for the solvent molecule probe radius c diff --git a/source/faces.f b/source/faces.f index 7b3c7bed9..1f6bb7a15 100644 --- a/source/faces.f +++ b/source/faces.f @@ -29,26 +29,40 @@ c c module faces + use sizes implicit none integer maxcls,maxtt - integer maxt,maxp - integer maxv,maxen - integer maxfn,maxc + integer maxt,maxp,maxv + integer maxen,maxfn,maxc integer maxep,maxfs integer maxcy,mxcyep integer maxfp,mxfpcy + parameter (maxcls=50*maxatm) + parameter (maxtt=25*maxatm) + parameter (maxt=3*maxatm) + parameter (maxp=2*maxatm) + parameter (maxv=5*maxatm) + parameter (maxen=5*maxatm) + parameter (maxfn=2*maxatm) + parameter (maxc=5*maxatm) + parameter (maxep=5*maxatm) + parameter (maxfs=3*maxatm) + parameter (maxcy=maxatm) + parameter (mxcyep=30) + parameter (maxfp=maxatm) + parameter (mxfpcy=10) c c c na number of atoms c pr probe radius c ar atomic radii -c axyz atomic coordinates +c a atomic coordinates c c integer na real*8 pr - real*8, allocatable :: ar(:) - real*8, allocatable :: axyz(:,:) + real*8 ar(maxatm) + real*8 a(3,maxatm) c c c skip if true, atom is not used @@ -57,10 +71,10 @@ module faces c abur atom buried c c - logical, allocatable :: skip(:) - logical, allocatable :: nosurf(:) - logical, allocatable :: afree(:) - logical, allocatable :: abur(:) + logical skip(maxatm) + logical nosurf(maxatm) + logical afree(maxatm) + logical abur(maxatm) c c c cls atom numbers of neighbors @@ -68,9 +82,9 @@ module faces c acls begin and end pointers for atoms neighbors c c - integer, allocatable :: cls(:) - integer, allocatable :: clst(:) - integer, allocatable :: acls(:,:) + integer cls(maxcls) + integer clst(maxcls) + integer acls(2,maxatm) c c c ntt number of temporary tori @@ -83,12 +97,12 @@ module faces c c integer ntt - integer, allocatable :: ttfe(:) - integer, allocatable :: ttle(:) - integer, allocatable :: enext(:) - integer, allocatable :: tta(:,:) - logical, allocatable :: ttbur(:) - logical, allocatable :: ttfree(:) + integer ttfe(maxtt) + integer ttle(maxtt) + integer enext(maxen) + integer tta(2,maxtt) + logical ttbur(maxtt) + logical ttfree(maxtt) c c c nt number of tori @@ -101,12 +115,12 @@ module faces c c integer nt - integer, allocatable :: tfe(:) - integer, allocatable :: ta(:,:) - real*8, allocatable :: tr(:) - real*8, allocatable :: t(:,:) - real*8, allocatable :: tax(:,:) - logical, allocatable :: tfree(:) + integer tfe(maxt) + integer ta(2,maxt) + real*8 tr(maxt) + real*8 t(3,maxt) + real*8 tax(3,maxt) + logical tfree(maxt) c c c np number of probe positions @@ -115,20 +129,20 @@ module faces c c integer np - integer, allocatable :: pa(:,:) - real*8, allocatable :: p(:,:) + integer pa(3,maxp) + real*8 p(3,maxp) c c +c v vertex coordinates c nv number of vertices c va vertex atom number c vp vertex probe number -c vxyz vertex coordinates c c integer nv - integer, allocatable :: va(:) - integer, allocatable :: vp(:) - real*8, allocatable :: vxyz(:,:) + integer va(maxv) + integer vp(maxv) + real*8 v(3,maxv) c c c nen number of concave edges @@ -139,8 +153,8 @@ module faces c integer nen integer nfn - integer, allocatable :: env(:,:) - integer, allocatable :: fnen(:,:) + integer env(2,maxen) + integer fnen(3,maxfn) c c c nc number of circles @@ -151,10 +165,10 @@ module faces c c integer nc - integer, allocatable :: ca(:) - integer, allocatable :: ct(:) - real*8, allocatable :: cr(:) - real*8, allocatable :: c(:,:) + integer ca(maxc) + integer ct(maxc) + real*8 cr(maxc) + real*8 c(3,maxc) c c c nep number of convex edges @@ -166,11 +180,11 @@ module faces c c integer nep - integer, allocatable :: epc(:) - integer, allocatable :: epv(:,:) - integer, allocatable :: afe(:) - integer, allocatable :: ale(:) - integer, allocatable :: epnext(:) + integer epc(maxep) + integer epv(2,maxep) + integer afe(maxatm) + integer ale(maxatm) + integer epnext(maxep) c c c nfs number of saddle faces @@ -179,8 +193,8 @@ module faces c c integer nfs - integer, allocatable :: fsen(:,:) - integer, allocatable :: fsep(:,:) + integer fsen(2,maxfs) + integer fsep(2,maxfs) c c c ncy number of cycles @@ -189,8 +203,8 @@ module faces c c integer ncy - integer, allocatable :: cynep(:) - integer, allocatable :: cyep(:,:) + integer cynep(maxcy) + integer cyep(mxcyep,maxcy) c c c nfp number of convex faces @@ -200,8 +214,8 @@ module faces c c integer nfp - integer, allocatable :: fpa(:) - integer, allocatable :: fpncy(:) - integer, allocatable :: fpcy(:,:) + integer fpa(maxfp) + integer fpncy(maxfp) + integer fpcy(mxfpcy,maxfp) save end diff --git a/source/final.f b/source/final.f index b27535d31..0e8c66811 100644 --- a/source/final.f +++ b/source/final.f @@ -33,8 +33,6 @@ subroutine final use deriv use dipole use disgeo - use domega - use faces use fracs use freeze use group @@ -47,7 +45,6 @@ subroutine final use light use merck use molcul - use moldyn use mpole use mutant use neigh @@ -66,7 +63,6 @@ subroutine final use qmstuf use refer use restrn - use rgddyn use rigid use ring use socket @@ -84,7 +80,6 @@ subroutine final use usolve use vdw use vibs - use warp implicit none c c @@ -215,82 +210,6 @@ subroutine final if (allocated(dbnd)) deallocate (dbnd) if (allocated(georad)) deallocate (georad) c -c deallocation of global arrays from module domega -c - if (allocated(tesum)) deallocate (tesum) - if (allocated(teb)) deallocate (teb) - if (allocated(tea)) deallocate (tea) - if (allocated(teba)) deallocate (teba) - if (allocated(teub)) deallocate (teub) - if (allocated(teaa)) deallocate (teaa) - if (allocated(teopb)) deallocate (teopb) - if (allocated(teopd)) deallocate (teopd) - if (allocated(teid)) deallocate (teid) - if (allocated(teit)) deallocate (teit) - if (allocated(tet)) deallocate (tet) - if (allocated(tept)) deallocate (tept) - if (allocated(tebt)) deallocate (tebt) - if (allocated(teat)) deallocate (teat) - if (allocated(tett)) deallocate (tett) - if (allocated(tev)) deallocate (tev) - if (allocated(tec)) deallocate (tec) - if (allocated(tecd)) deallocate (tecd) - if (allocated(ted)) deallocate (ted) - if (allocated(tem)) deallocate (tem) - if (allocated(tep)) deallocate (tep) - if (allocated(ter)) deallocate (ter) - if (allocated(tes)) deallocate (tes) - if (allocated(telf)) deallocate (telf) - if (allocated(teg)) deallocate (teg) - if (allocated(tex)) deallocate (tex) -c -c deallocation of global arrays from module faces -c - if (allocated(ar)) deallocate (ar) - if (allocated(axyz)) deallocate (axyz) - if (allocated(skip)) deallocate (skip) - if (allocated(nosurf)) deallocate (nosurf) - if (allocated(afree)) deallocate (afree) - if (allocated(abur)) deallocate (abur) - if (allocated(cls)) deallocate (cls) - if (allocated(clst)) deallocate (clst) - if (allocated(acls)) deallocate (acls) - if (allocated(ttfe)) deallocate (ttfe) - if (allocated(ttle)) deallocate (ttle) - if (allocated(enext)) deallocate (enext) - if (allocated(tta)) deallocate (tta) - if (allocated(ttbur)) deallocate (ttbur) - if (allocated(ttfree)) deallocate (ttfree) - if (allocated(tfe)) deallocate (tfe) - if (allocated(ta)) deallocate (ta) - if (allocated(tr)) deallocate (tr) - if (allocated(t)) deallocate (t) - if (allocated(tax)) deallocate (tax) - if (allocated(tfree)) deallocate (tfree) - if (allocated(pa)) deallocate (pa) - if (allocated(p)) deallocate (p) - if (allocated(va)) deallocate (va) - if (allocated(vp)) deallocate (vp) - if (allocated(vxyz)) deallocate (vxyz) - if (allocated(env)) deallocate (env) - if (allocated(fnen)) deallocate (fnen) - if (allocated(ca)) deallocate (ca) - if (allocated(ct)) deallocate (ct) - if (allocated(cr)) deallocate (cr) - if (allocated(c)) deallocate (c) - if (allocated(epc)) deallocate (epc) - if (allocated(epv)) deallocate (epv) - if (allocated(afe)) deallocate (afe) - if (allocated(ale)) deallocate (ale) - if (allocated(epnext)) deallocate (epnext) - if (allocated(fsen)) deallocate (fsen) - if (allocated(fsep)) deallocate (fsep) - if (allocated(cynep)) deallocate (cynep) - if (allocated(cyep)) deallocate (cyep) - if (allocated(fpa)) deallocate (fpa) - if (allocated(fpncy)) deallocate (fpncy) - if (allocated(fpcy)) deallocate (fpcy) -c c deallocation of global arrays from module fracs c if (allocated(xfrac)) deallocate (xfrac) @@ -405,12 +324,6 @@ subroutine final if (allocated(molcule)) deallocate (molcule) if (allocated(molmass)) deallocate (molmass) c -c deallocation of global arrays from module moldyn -c - if (allocated(v)) deallocate (v) - if (allocated(a)) deallocate (a) - if (allocated(aalt)) deallocate (aalt) -c c deallocation of global arrays from module mpole c if (allocated(ipole)) deallocate (ipole) @@ -584,18 +497,6 @@ subroutine final if (allocated(gfix)) deallocate (gfix) if (allocated(chir)) deallocate (chir) c -c deallocation of global arrays from module rgddyn -c - if (allocated(xcmo)) deallocate (xcmo) - if (allocated(ycmo)) deallocate (ycmo) - if (allocated(zcmo)) deallocate (zcmo) - if (allocated(vcm)) deallocate (vcm) - if (allocated(wcm)) deallocate (wcm) - if (allocated(lm)) deallocate (lm) - if (allocated(vc)) deallocate (vc) - if (allocated(wc)) deallocate (wc) - if (allocated(linear)) deallocate (linear) -c c deallocation of global arrays from module rigid c if (allocated(xrb)) deallocate (xrb) @@ -610,32 +511,14 @@ subroutine final if (allocated(iring5)) deallocate (iring5) if (allocated(iring6)) deallocate (iring6) c -c deallocation of global arrays from module solute -c - if (allocated(rsolv)) deallocate (rsolv) - if (allocated(asolv)) deallocate (asolv) - if (allocated(rborn)) deallocate (rborn) - if (allocated(drb)) deallocate (drb) - if (allocated(drbp)) deallocate (drbp) - if (allocated(drobc)) deallocate (drobc) - if (allocated(gpol)) deallocate (gpol) - if (allocated(shct)) deallocate (shct) - if (allocated(aobc)) deallocate (aobc) - if (allocated(bobc)) deallocate (bobc) - if (allocated(gobc)) deallocate (gobc) - if (allocated(vsolv)) deallocate (vsolv) - if (allocated(wace)) deallocate (wace) - if (allocated(s2ace)) deallocate (s2ace) - if (allocated(uace)) deallocate (uace) -c c deallocation of global arrays from module stodyn c if (allocated(fgamma)) deallocate (fgamma) c c deallocation of global arrays from module strbnd c - if (allocated(isb)) deallocate (isb) - if (allocated(sbk)) deallocate (sbk) + if (allocated(ist)) deallocate (isb) + if (allocated(kst)) deallocate (sbk) c c deallocation of global arrays from module strtor c @@ -709,10 +592,6 @@ subroutine final if (allocated(phik)) deallocate (phik) if (allocated(pwork)) deallocate (pwork) c -c deallocation of global arrays from module warp -c - if (allocated(m2)) deallocate (m2) -c c free memory used by the APBS Poisson-Boltzmann solver c if (solvtyp .eq. 'PB') then diff --git a/source/flatten.f b/source/flatten.f index 07787c9fc..fae95d1a1 100644 --- a/source/flatten.f +++ b/source/flatten.f @@ -112,12 +112,6 @@ subroutine flatten end if end if c -c perform dynamic allocation of some global arrays -c - if (use_gda) then - if (.not. allocated(m2)) allocate (m2(n)) - end if -c c set second moment of Gaussian on each atom for GDA methods c if (use_gda) then diff --git a/source/getkey.f b/source/getkey.f index d04cc0ca9..0c22e7407 100644 --- a/source/getkey.f +++ b/source/getkey.f @@ -17,6 +17,7 @@ c c subroutine getkey + use sizes use argue use files use iounit diff --git a/source/getprm.f b/source/getprm.f index 1627a7b94..b686ae989 100644 --- a/source/getprm.f +++ b/source/getprm.f @@ -17,6 +17,7 @@ c c subroutine getprm + use sizes use files use iounit use keys diff --git a/source/gradient.f b/source/gradient.f index 4dbda2a96..c391a3060 100644 --- a/source/gradient.f +++ b/source/gradient.f @@ -70,32 +70,32 @@ subroutine gradient (energy,derivs) c if (first) then first = .false. - if (.not. allocated(desum)) allocate (desum(3,n)) - if (.not. allocated(deb)) allocate (deb(3,n)) - if (.not. allocated(dea)) allocate (dea(3,n)) - if (.not. allocated(deba)) allocate (deba(3,n)) - if (.not. allocated(deub)) allocate (deub(3,n)) - if (.not. allocated(deaa)) allocate (deaa(3,n)) - if (.not. allocated(deopb)) allocate (deopb(3,n)) - if (.not. allocated(deopd)) allocate (deopd(3,n)) - if (.not. allocated(deid)) allocate (deid(3,n)) - if (.not. allocated(deit)) allocate (deit(3,n)) - if (.not. allocated(det)) allocate (det(3,n)) - if (.not. allocated(dept)) allocate (dept(3,n)) - if (.not. allocated(debt)) allocate (debt(3,n)) - if (.not. allocated(deat)) allocate (deat(3,n)) - if (.not. allocated(dett)) allocate (dett(3,n)) - if (.not. allocated(dev)) allocate (dev(3,n)) - if (.not. allocated(dec)) allocate (dec(3,n)) - if (.not. allocated(decd)) allocate (decd(3,n)) - if (.not. allocated(ded)) allocate (ded(3,n)) - if (.not. allocated(dem)) allocate (dem(3,n)) - if (.not. allocated(dep)) allocate (dep(3,n)) - if (.not. allocated(der)) allocate (der(3,n)) - if (.not. allocated(des)) allocate (des(3,n)) - if (.not. allocated(delf)) allocate (delf(3,n)) - if (.not. allocated(deg)) allocate (deg(3,n)) - if (.not. allocated(dex)) allocate (dex(3,n)) + if (.not. allocated(desum)) allocate (desum(3,maxatm)) + if (.not. allocated(deb)) allocate (deb(3,maxatm)) + if (.not. allocated(dea)) allocate (dea(3,maxatm)) + if (.not. allocated(deba)) allocate (deba(3,maxatm)) + if (.not. allocated(deub)) allocate (deub(3,maxatm)) + if (.not. allocated(deaa)) allocate (deaa(3,maxatm)) + if (.not. allocated(deopb)) allocate (deopb(3,maxatm)) + if (.not. allocated(deopd)) allocate (deopd(3,maxatm)) + if (.not. allocated(deid)) allocate (deid(3,maxatm)) + if (.not. allocated(deit)) allocate (deit(3,maxatm)) + if (.not. allocated(det)) allocate (det(3,maxatm)) + if (.not. allocated(dept)) allocate (dept(3,maxatm)) + if (.not. allocated(debt)) allocate (debt(3,maxatm)) + if (.not. allocated(deat)) allocate (deat(3,maxatm)) + if (.not. allocated(dett)) allocate (dett(3,maxatm)) + if (.not. allocated(dev)) allocate (dev(3,maxatm)) + if (.not. allocated(dec)) allocate (dec(3,maxatm)) + if (.not. allocated(decd)) allocate (decd(3,maxatm)) + if (.not. allocated(ded)) allocate (ded(3,maxatm)) + if (.not. allocated(dem)) allocate (dem(3,maxatm)) + if (.not. allocated(dep)) allocate (dep(3,maxatm)) + if (.not. allocated(der)) allocate (der(3,maxatm)) + if (.not. allocated(des)) allocate (des(3,maxatm)) + if (.not. allocated(delf)) allocate (delf(3,maxatm)) + if (.not. allocated(deg)) allocate (deg(3,maxatm)) + if (.not. allocated(dex)) allocate (dex(3,maxatm)) end if c c zero out each of the first derivative components diff --git a/source/gradrot.f b/source/gradrot.f index 6594c6eda..dcf3f94ac 100644 --- a/source/gradrot.f +++ b/source/gradrot.f @@ -33,42 +33,7 @@ subroutine gradrot (energy,derivs) real*8 xterm,yterm,zterm real*8 derivs(*) real*8, allocatable :: g(:,:) - logical first - save first - data first / .true. / -c -c -c perform dynamic allocation of some global arrays -c - if (first) then - first = .false. - if (.not. allocated(tesum)) allocate (tesum(nomega)) - if (.not. allocated(teb)) allocate (teb(nomega)) - if (.not. allocated(tea)) allocate (tea(nomega)) - if (.not. allocated(teba)) allocate (teba(nomega)) - if (.not. allocated(teub)) allocate (teub(nomega)) - if (.not. allocated(teaa)) allocate (teaa(nomega)) - if (.not. allocated(teopb)) allocate (teopb(nomega)) - if (.not. allocated(teopd)) allocate (teopd(nomega)) - if (.not. allocated(teid)) allocate (teid(nomega)) - if (.not. allocated(teit)) allocate (teit(nomega)) - if (.not. allocated(tet)) allocate (tet(nomega)) - if (.not. allocated(tept)) allocate (tept(nomega)) - if (.not. allocated(tebt)) allocate (tebt(nomega)) - if (.not. allocated(teat)) allocate (teat(nomega)) - if (.not. allocated(tett)) allocate (tett(nomega)) - if (.not. allocated(tev)) allocate (tev(nomega)) - if (.not. allocated(tec)) allocate (tec(nomega)) - if (.not. allocated(tecd)) allocate (tecd(nomega)) - if (.not. allocated(ted)) allocate (ted(nomega)) - if (.not. allocated(tem)) allocate (tem(nomega)) - if (.not. allocated(tep)) allocate (tep(nomega)) - if (.not. allocated(ter)) allocate (ter(nomega)) - if (.not. allocated(tes)) allocate (tes(nomega)) - if (.not. allocated(telf)) allocate (telf(nomega)) - if (.not. allocated(teg)) allocate (teg(nomega)) - if (.not. allocated(tex)) allocate (tex(nomega)) - end if +c c c zero out individual components of torsional derivatives c diff --git a/source/hessian.f b/source/hessian.f index 1f5ab4814..c14870336 100644 --- a/source/hessian.f +++ b/source/hessian.f @@ -109,9 +109,9 @@ subroutine hessian (h,hinit,hstop,hindex,hdiag) c if (first) then first = .false. - if (.not. allocated(hessx)) allocate (hessx(3,n)) - if (.not. allocated(hessy)) allocate (hessy(3,n)) - if (.not. allocated(hessz)) allocate (hessz(3,n)) + if (.not. allocated(hessx)) allocate (hessx(3,maxatm)) + if (.not. allocated(hessy)) allocate (hessy(3,maxatm)) + if (.not. allocated(hessz)) allocate (hessz(3,maxatm)) end if c c zero out the Hessian elements for the current atom diff --git a/source/initatom.f b/source/initatom.f index 66b91cce3..f3e3fb3cc 100644 --- a/source/initatom.f +++ b/source/initatom.f @@ -17,6 +17,7 @@ c c subroutine initatom + use sizes use ptable implicit none integer i diff --git a/source/initprm.f b/source/initprm.f index 70ed36e14..046017f82 100644 --- a/source/initprm.f +++ b/source/initprm.f @@ -55,7 +55,7 @@ subroutine initprm use units use vdwpot implicit none - integer i,j + integer i,j,k character*3 blank3 character*8 blank8 character*12 blank12 @@ -296,7 +296,7 @@ subroutine initprm c poltyp = 'MUTUAL' politer = 500 - poleps = 0.00001d0 + poleps = 0.000001d0 udiag = 2.0d0 d1scale = 0.0d0 d2scale = 1.0d0 diff --git a/source/initres.f b/source/initres.f index ff76b7794..1b343b09e 100644 --- a/source/initres.f +++ b/source/initres.f @@ -17,6 +17,7 @@ c c subroutine initres + use sizes use resdue implicit none integer i diff --git a/source/isokinetic.f b/source/isokinetic.f new file mode 100644 index 000000000..cd081f610 --- /dev/null +++ b/source/isokinetic.f @@ -0,0 +1,476 @@ +c +c +c ################################################### +c ## COPYRIGHT (C) 2011 by Jay William Ponder ## +c ## All Rights Reserved ## +c ################################################### +c +c ############################################################# +c ## ## +c ## subroutine isokinetic ## +c ## ## +c ############################################################# +c +c +c + subroutine isokinetic (istep,dt) + use sizes + use atomid + use freeze + use moldyn + use units + use virial + use atoms + use bath + use bound + use inform + use iounit + use keys + use mdstuf + use potent + use solute + use stodyn + use usage + implicit none + integer i,j,k + integer istep + integer nalt + real*8 dt,dt_2,tempstor,temp2 + real*8 dta,dta_2,kBT,LkT + real*8 epot,etot + real*8 eksum,eps + real*8 temp,pres + real*8 ealt,dalt + real*8 ekin(3,3) + real*8 stress(3,3) + real*8 viralt(3,3) + real*8, allocatable :: xold(:) + real*8, allocatable :: yold(:) + real*8, allocatable :: zold(:) + real*8, allocatable :: derivs(:,:) + +c +c set some values for the dynamics integration +c + dt_2 = 0.5d0 * dt + dta = dt / (dble(nrespa)) + dta_2 = 0.5d0 * dta + + kBT = boltzmann*kelvin + LkT = len_nhc*kBT +c +c perform dynamic allocation of some local arrays +c + allocate (xold(n)) + allocate (yold(n)) + allocate (zold(n)) + allocate (derivs(3,n)) + +c +c initialize virial from fast-evolving potential energy terms +c + do i = 1, 3 + do j = 1, 3 + viralt(j,i) = 0.0d0 + end do + end do + +c +c +c begin inner RESPA loop +c + do k = 1, nrespa + do i = 1, n + if (use(i)) then + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + end if + end do + +c +c +c nose-like isokinetic operator +c + call appisok (dta) +c +c +c force-like isokinetic operator +c + + call appisok2 (dta) +c +c Evolve positions +c + if (use_rattle) then + do i = 1, n + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end do + call rattle (dta,xold,yold,zold) + else + do i = 1, n + if (use(i)) then + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end if + end do + end if +c +c get the potential energy and atomic forces +c + if (k.eq.nrespa) then + call gradfast (ealt,derivs) + + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + end if + + else + call gradfast (ealt,derivs) + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + end if + end if + +c +c +c force-like isokinetic operator +c + call appisok2 (dta) +c +c +c nose-like isokinetic operator +c + call appisok (dta) + + end do + +c +c find the constraint-corrected full-step velocities +c + if (use_rattle) call rattle2 (dt) +c +c total potential and virial from sum of fast and slow parts +c + epot = epot + ealt + do i = 1, 3 + do j = 1, 3 + vir(j,i) = vir(j,i) + viralt(j,i) + end do + end do + + +c +c perform deallocation of some local arrays +c + + deallocate (derivs) + deallocate (xold) + deallocate (yold) + deallocate (zold) +c +c total energy is sum of kinetic and potential energies +c + etot = eksum + epot + + call kinetic (eksum,ekin) + call pressure (dt,epot,ekin,temp,pres,stress) + temp = 2.0d0 * eksum / (dble(nfree) * gasconst) + + +c tempstor=0.0d0 +c do i = 1, n +c tempstor=tempstor+(mass(i)*v(1,i)*v(1,i)) +c temp2=0.5*q_iso1(1,1,i)*v_iso1(1,1,i)*v_iso1(1,1,i) +c tempstor=tempstor+temp2 +c +c tempstor=tempstor+(mass(i)*v(2,i)*v(2,i)) +c temp2=0.5*q_iso1(1,2,i)*v_iso1(1,2,i)*v_iso1(1,2,i) +c tempstor=tempstor+temp2 +c +c tempstor=tempstor+(mass(i)*v(3,i)*v(3,i)) +c temp2=0.5*q_iso1(1,3,i)*v_iso1(1,3,i)*v_iso1(1,3,i) +c tempstor=tempstor+temp2 +c end do +c write (iout,10) tempstor/(3.0*n*kBT) +c 10 format (3f16.14) + + +c +c compute statistics and save trajectory for this step +c + call mdstat (istep,dt,etot,epot,eksum,temp,pres) + call mdsave (istep,dt,epot) + call mdrest (istep) + + + + + return + end + + +c +c ################################################################## +c ## ## +c ## subroutine applyisok ## +c ## ## +c ################################################################## +c +c +c +c +c + subroutine appisok (dt) + use sizes + use atomid + use freeze + use moldyn + use units + use virial + use atoms + use bath + use bound + use inform + use iounit + use keys + use mdstuf + use potent + use solute + use stodyn + use usage + implicit none + integer nc,ns,i,j,k,iatom,inhc,ixyz,index + integer len_m1,len_p1,len_m2 + real*8 dt,dtc,dts + real*8 dt2,dt4,dt8 + real*8 tempstor,kBT,LkT,temp2,temp3 + real*8 w(3) + + len_m1=len_nhc-1 + len_p1=len_nhc+1 + len_m2=len_nhc-2 + + kBT = boltzmann*kelvin + LkT = len_nhc*kBT + + nc = 5 + ns = 3 + dtc = dt / dble(nc) + w(1) = 1.0d0 / (2.0d0-2.0d0**(1.0d0/3.0d0)) + w(2) = 1.0d0 - 2.0d0*w(1) + w(3) = w(1) + + +c +c use multiple time steps and Suzuki-Yoshida decomposition +c + do k = 1, nc + do j = 1, ns + + dts = w(j) * dtc + dt2 = 0.5d0 * dts + dt4 = 0.25d0 * dts + dt8 = 0.125d0 * dts + + do iatom = 1, n + do ixyz = 1, 3 + do inhc = 2, len_nhc + f_iso1(inhc,ixyz,iatom)=((q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom))-kBT) + end do + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + v_iso1(len_nhc,ixyz,iatom)=v_iso1(len_nhc,ixyz,iatom) + & +dt4*f_iso1(len_nhc,ixyz,iatom) + & /q_iso1(len_nhc,ixyz,iatom) + end do + end do + + + do iatom = 1, n + do inhc = len_m1, 2, -1 + do ixyz = 1, 3 + index=inhc+1 + tempstor = exp(-dt8*v_iso1(index,ixyz,iatom)) + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & *tempstor*tempstor + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & +dt4*tempstor*f_iso1(inhc,ixyz,iatom) + & /q_iso1(index,ixyz,iatom) + end do + end do + end do + + + do iatom = 1, n + do ixyz = 1, 3 + temp2=v(ixyz,iatom)*v(ixyz,iatom)*mass(iatom) + & +(0.5*q_iso1(1,ixyz,iatom)*v_iso1(1,ixyz,iatom) + & *v_iso1(1,ixyz,iatom) + & *exp(-2.0*dt2*v_iso1(2,ixyz,iatom))) + temp3=sqrt(kBT/temp2) + v(ixyz,iatom)=v(ixyz,iatom)*temp3 + v_iso1(1,ixyz,iatom)=v_iso1(1,ixyz,iatom) + & *temp3*exp(-dt2*v_iso1(2,ixyz,iatom)) + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + do inhc = 2, len_nhc + f_iso1(inhc,ixyz,iatom)=((q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom) + & *v_iso1(inhc-1,ixyz,iatom))-kBT) + end do + end do + end do + + do iatom = 1, n + do inhc = 2, len_m1 + do ixyz = 1, 3 + index=inhc+1 + tempstor = exp(-dt8*v_iso1(index,ixyz,iatom)) + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & *tempstor*tempstor + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) + & +dt4*tempstor*f_iso1(inhc,ixyz,iatom) + & /q_iso1(index,ixyz,iatom) + end do + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + v_iso1(len_nhc,ixyz,iatom)=v_iso1(len_nhc,ixyz,iatom) + & +dt4*f_iso1(len_nhc,ixyz,iatom) + & /q_iso1(len_nhc,ixyz,iatom) + end do + end do + + end do + end do + + + return + end +c +c ################################################################## +c ## ## +c ## subroutine applyisok2 ## +c ## ## +c ################################################################## +c +c +c + subroutine appisok2 (dt) + use sizes + use atomid + use freeze + use moldyn + use units + use virial + use atoms + use bath + use bound + use inform + use iounit + use keys + use mdstuf + use potent + use solute + use stodyn + use usage + implicit none + real*8 kBT,LkT,isoA,isoB + real*8 dt,dt_2,rootB,arg + real*8 iso_s,iso_sdot + real*8 epot,etot + integer iatom,inhc,ixyz,i,j,k + + dt_2 = 0.5d0*dt + kBT = boltzmann*kelvin + LkT = len_nhc*kBT + + + do iatom = 1, n + do ixyz = 1, 3 + isoA = mass(iatom)*a(ixyz,iatom)*v(ixyz,iatom)/kBT + isoB = mass(iatom)*a(ixyz,iatom)*a(ixyz,iatom)/kBT + rootB = sqrt(isoB) + arg = dt_2*rootB + + if (arg > 0.00001) then + iso_s = (1.0/rootB)*sinh(arg) + + & (isoA/isoB)*(cosh(arg)-1.0) + iso_sdot = cosh(arg) + (isoA/rootB)*sinh(arg) + else + iso_s = ((((isoB*isoA/24.0)*dt_2 +isoB/6.0)*dt_2 + + & 0.5*isoA)*dt_2 + 1.0)*dt_2 + iso_sdot = (((isoB*isoA/6.0)*dt_2 +0.5*isoB)*dt_2 + + & isoA)*dt_2 + 1.0 + end if + + v(ixyz,iatom) = (v(ixyz,iatom) + ((a(ixyz,iatom))*iso_s)) + & /iso_sdot + v_iso1(1,ixyz,iatom)=v_iso1(1,ixyz,iatom)/iso_sdot + end do + end do + + + return + end diff --git a/source/kangang.f b/source/kangang.f index a5ce88ff0..a66ffe335 100644 --- a/source/kangang.f +++ b/source/kangang.f @@ -29,12 +29,10 @@ subroutine kangang use kanang use keys use potent - use tors implicit none integer i,j,k,m,next integer it,ia,ic integer nang,jang,kang - integer maxaa real*8 fa,faa,aak(3) logical header character*20 keyword @@ -77,9 +75,8 @@ subroutine kangang c c perform dynamic allocation of some global arrays c - maxaa = 15 * n - if (.not. allocated(iaa)) allocate (iaa(2,maxaa)) - if (.not. allocated(kaa)) allocate (kaa(maxaa)) + if (.not. allocated(iaa)) allocate (iaa(2,maxtors)) + if (.not. allocated(kaa)) allocate (kaa(maxtors)) c c assign the angle-angle parameters for each angle pair c diff --git a/source/kangle.f b/source/kangle.f index 9c95eca70..c1aa9217e 100644 --- a/source/kangle.f +++ b/source/kangle.f @@ -307,10 +307,10 @@ subroutine kangle c c perform dynamic allocation of some global arrays c - if (.not. allocated(ak)) allocate (ak(nangle)) - if (.not. allocated(anat)) allocate (anat(nangle)) - if (.not. allocated(afld)) allocate (afld(nangle)) - if (.not. allocated(angtyp)) allocate (angtyp(nangle)) + if (.not. allocated(ak)) allocate (ak(maxang)) + if (.not. allocated(anat)) allocate (anat(maxang)) + if (.not. allocated(afld)) allocate (afld(maxang)) + if (.not. allocated(angtyp)) allocate (angtyp(maxang)) c c use special angle parameter assignment method for MMFF c @@ -685,7 +685,7 @@ subroutine kanglem if (c(ina) .eq. 1000.0d0) goto 20 if (c(inb) .eq. 1000.0d0) goto 20 if (c(inc) .eq. 1000.0d0) goto 20 - do k = 1, nbond + do k = 1, maxbnd if ((min(ia,ib).eq.ibnd(1,k)) .and. & (max(ia,ib).eq.ibnd(2,k))) then bnd_ab = k @@ -739,7 +739,7 @@ subroutine kanglem end if if (ring3) anat(i) = 60.0d0 if (ring4) anat(i) = 90.0d0 - do k = 1, nbond + do k = 1, maxbnd if ((min(ia,ib).eq.ibnd(1,k)) .and. & (max(ia,ib).eq.ibnd(2,k))) then bnd_ab = k diff --git a/source/kangtor.f b/source/kangtor.f index 169891dee..579d8168c 100644 --- a/source/kangtor.f +++ b/source/kangtor.f @@ -140,8 +140,8 @@ subroutine kangtor c c perform dynamic allocation of some global arrays c - if (.not. allocated(iat)) allocate (iat(3,ntors)) - if (.not. allocated(kant)) allocate (kant(6,ntors)) + if (.not. allocated(iat)) allocate (iat(3,maxtors)) + if (.not. allocated(kant)) allocate (kant(6,maxtors)) c c assign the angle-torsion parameters for each torsion c diff --git a/source/kbond.f b/source/kbond.f index c994f3d4d..5f70555d9 100644 --- a/source/kbond.f +++ b/source/kbond.f @@ -178,8 +178,8 @@ subroutine kbond c c perform dynamic allocation of some global arrays c - if (.not. allocated(bk)) allocate (bk(nbond)) - if (.not. allocated(bl)) allocate (bl(nbond)) + if (.not. allocated(bk)) allocate (bk(maxbnd)) + if (.not. allocated(bl)) allocate (bl(maxbnd)) c c use special bond parameter assignment method for MMFF c diff --git a/source/kcharge.f b/source/kcharge.f index e962c478e..5115031c0 100644 --- a/source/kcharge.f +++ b/source/kcharge.f @@ -80,11 +80,11 @@ subroutine kcharge c c perform dynamic allocation of some global arrays c - if (.not. allocated(iion)) allocate (iion(n)) - if (.not. allocated(jion)) allocate (jion(n)) - if (.not. allocated(kion)) allocate (kion(n)) - if (.not. allocated(chglist)) allocate (chglist(n)) - if (.not. allocated(pchg)) allocate (pchg(n)) + if (.not. allocated(iion)) allocate (iion(maxatm)) + if (.not. allocated(jion)) allocate (jion(maxatm)) + if (.not. allocated(kion)) allocate (kion(maxatm)) + if (.not. allocated(chglist)) allocate (chglist(maxatm)) + if (.not. allocated(pchg)) allocate (pchg(maxatm)) c c find and store all the atomic partial charges c diff --git a/source/kdipole.f b/source/kdipole.f index d397714ab..ce9d36410 100644 --- a/source/kdipole.f +++ b/source/kdipole.f @@ -196,9 +196,9 @@ subroutine kdipole c c perform dynamic allocation of some global arrays c - if (.not. allocated(idpl)) allocate (idpl(2,nbond)) - if (.not. allocated(bdpl)) allocate (bdpl(nbond)) - if (.not. allocated(sdpl)) allocate (sdpl(nbond)) + if (.not. allocated(idpl)) allocate (idpl(2,maxbnd)) + if (.not. allocated(bdpl)) allocate (bdpl(maxbnd)) + if (.not. allocated(sdpl)) allocate (sdpl(maxbnd)) c c find and store all the bond dipole moments c diff --git a/source/keys.f b/source/keys.f index 1a4e23e3c..f981a807c 100644 --- a/source/keys.f +++ b/source/keys.f @@ -12,16 +12,13 @@ c ############################################################# c c -c maxkey maximum number of lines in the keyword file -c c nkey number of nonblank lines in the keyword file c keyline contents of each individual keyword file line c c module keys + use sizes implicit none - integer maxkey - parameter (maxkey=25000) integer nkey character*120 keyline(maxkey) save diff --git a/source/kimprop.f b/source/kimprop.f index 9f0e2b834..f58aa8610 100644 --- a/source/kimprop.f +++ b/source/kimprop.f @@ -27,7 +27,6 @@ subroutine kimprop use keys use kiprop use potent - use tors implicit none integer i,j,k,ndi integer ia,ib,ic,id @@ -109,9 +108,9 @@ subroutine kimprop c c perform dynamic allocation of some global arrays c - if (.not. allocated(iiprop)) allocate (iiprop(4,ntors)) - if (.not. allocated(kprop)) allocate (kprop(ntors)) - if (.not. allocated(vprop)) allocate (vprop(ntors)) + if (.not. allocated(iiprop)) allocate (iiprop(4,maxtors)) + if (.not. allocated(kprop)) allocate (kprop(maxtors)) + if (.not. allocated(vprop)) allocate (vprop(maxtors)) c c assign improper dihedral parameters for each improper angle; c multiple symmetrical parameters are given partial weights diff --git a/source/kimptor.f b/source/kimptor.f index da060f3f7..7629931ab 100644 --- a/source/kimptor.f +++ b/source/kimptor.f @@ -28,7 +28,6 @@ subroutine kimptor use kitors use math use potent - use tors implicit none integer i,j,k,nti integer ia,ib,ic,id @@ -120,10 +119,10 @@ subroutine kimptor c c perform dynamic allocation of some global arrays c - if (.not. allocated(iitors)) allocate (iitors(4,ntors)) - if (.not. allocated(itors1)) allocate (itors1(4,ntors)) - if (.not. allocated(itors2)) allocate (itors2(4,ntors)) - if (.not. allocated(itors3)) allocate (itors3(4,ntors)) + if (.not. allocated(iitors)) allocate (iitors(4,maxtors)) + if (.not. allocated(itors1)) allocate (itors1(4,maxtors)) + if (.not. allocated(itors2)) allocate (itors2(4,maxtors)) + if (.not. allocated(itors3)) allocate (itors3(4,maxtors)) c c assign improper torsional parameters for each improper torsion; c multiple symmetrical parameters are given partial weights diff --git a/source/kmpole.f b/source/kmpole.f index b4c6fcdc6..694d7abf8 100644 --- a/source/kmpole.f +++ b/source/kmpole.f @@ -206,19 +206,19 @@ subroutine kmpole c c perform dynamic allocation of some global arrays c - if (.not. allocated(ipole)) allocate (ipole(n)) - if (.not. allocated(polsiz)) allocate (polsiz(n)) - if (.not. allocated(pollist)) allocate (pollist(n)) - if (.not. allocated(zaxis)) allocate (zaxis(n)) - if (.not. allocated(xaxis)) allocate (xaxis(n)) - if (.not. allocated(yaxis)) allocate (yaxis(n)) - if (.not. allocated(pole)) allocate (pole(maxpole,n)) - if (.not. allocated(rpole)) allocate (rpole(maxpole,n)) - if (.not. allocated(polaxe)) allocate (polaxe(n)) - if (.not. allocated(np11)) allocate (np11(n)) - if (.not. allocated(np12)) allocate (np12(n)) - if (.not. allocated(np13)) allocate (np13(n)) - if (.not. allocated(np14)) allocate (np14(n)) + if (.not. allocated(ipole)) allocate (ipole(maxatm)) + if (.not. allocated(polsiz)) allocate (polsiz(maxatm)) + if (.not. allocated(pollist)) allocate (pollist(maxatm)) + if (.not. allocated(zaxis)) allocate (zaxis(maxatm)) + if (.not. allocated(xaxis)) allocate (xaxis(maxatm)) + if (.not. allocated(yaxis)) allocate (yaxis(maxatm)) + if (.not. allocated(pole)) allocate (pole(maxpole,maxatm)) + if (.not. allocated(rpole)) allocate (rpole(maxpole,maxatm)) + if (.not. allocated(polaxe)) allocate (polaxe(maxatm)) + if (.not. allocated(np11)) allocate (np11(maxatm)) + if (.not. allocated(np12)) allocate (np12(maxatm)) + if (.not. allocated(np13)) allocate (np13(maxatm)) + if (.not. allocated(np14)) allocate (np14(maxatm)) c c zero out local axes, multipoles and polarization attachments c @@ -514,16 +514,9 @@ subroutine kmpole polsiz(i) = size end do c -c perform dynamic allocation of some global arrays -c - if (.not. use_polar) then - if (.not. allocated(uind)) allocate (uind(3,n)) - if (.not. allocated(uinp)) allocate (uinp(3,n)) - if (.not. allocated(uinds)) allocate (uinds(3,n)) - if (.not. allocated(uinps)) allocate (uinps(3,n)) -c c if polarization not used, zero out induced dipoles c + if (.not. use_polar) then do i = 1, n do j = 1, 3 uind(j,i) = 0.0d0 diff --git a/source/kopbend.f b/source/kopbend.f index 305d3614e..7adead9f1 100644 --- a/source/kopbend.f +++ b/source/kopbend.f @@ -140,9 +140,9 @@ subroutine kopbend c c perform dynamic allocation of some global arrays c - if (.not. allocated(iopb)) allocate (iopb(nangle)) - if (.not. allocated(opbk)) allocate (opbk(nangle)) - if (.not. allocated(angtyp)) allocate (angtyp(nangle)) + if (.not. allocated(iopb)) allocate (iopb(maxang)) + if (.not. allocated(opbk)) allocate (opbk(maxang)) + if (.not. allocated(angtyp)) allocate (angtyp(maxang)) c c assign out-of-plane bending parameters for each angle c diff --git a/source/kopdist.f b/source/kopdist.f index d285ca163..fb2f5b48f 100644 --- a/source/kopdist.f +++ b/source/kopdist.f @@ -127,9 +127,9 @@ subroutine kopdist c c perform dynamic allocation of some global arrays c - if (.not. allocated(iopd)) allocate (iopd(4,n)) - if (.not. allocated(opdk)) allocate (opdk(n)) - if (.not. allocated(angtyp)) allocate (angtyp(nangle)) + if (.not. allocated(iopd)) allocate (iopd(4,maxatm)) + if (.not. allocated(opdk)) allocate (opdk(maxatm)) + if (.not. allocated(angtyp)) allocate (angtyp(maxang)) c c assign out-of-plane distance parameters for trigonal sites c diff --git a/source/kpitors.f b/source/kpitors.f index b09893256..91e94aff7 100644 --- a/source/kpitors.f +++ b/source/kpitors.f @@ -28,7 +28,6 @@ subroutine kpitors use kpitor use pitors use potent - use tors implicit none integer i,j,npt integer ia,ib @@ -102,8 +101,8 @@ subroutine kpitors c c perform dynamic allocation of some global arrays c - if (.not. allocated(ipit)) allocate (ipit(6,ntors)) - if (.not. allocated(kpit)) allocate (kpit(ntors)) + if (.not. allocated(ipit)) allocate (ipit(6,maxtors)) + if (.not. allocated(kpit)) allocate (kpit(maxtors)) c c assign pi-orbital torsion parameters as required c diff --git a/source/kpolar.f b/source/kpolar.f index d2ed98261..d164e843f 100644 --- a/source/kpolar.f +++ b/source/kpolar.f @@ -95,13 +95,13 @@ subroutine kpolar c c perform dynamic allocation of some global arrays c - if (.not. allocated(polarity)) allocate (polarity(n)) - if (.not. allocated(thole)) allocate (thole(n)) - if (.not. allocated(pdamp)) allocate (pdamp(n)) - if (.not. allocated(uind)) allocate (uind(3,n)) - if (.not. allocated(uinp)) allocate (uinp(3,n)) - if (.not. allocated(uinds)) allocate (uinds(3,n)) - if (.not. allocated(uinps)) allocate (uinps(3,n)) + if (.not. allocated(polarity)) allocate (polarity(maxatm)) + if (.not. allocated(thole)) allocate (thole(maxatm)) + if (.not. allocated(pdamp)) allocate (pdamp(maxatm)) + if (.not. allocated(uind)) allocate (uind(3,maxatm)) + if (.not. allocated(uinp)) allocate (uinp(3,maxatm)) + if (.not. allocated(uinds)) allocate (uinds(3,maxatm)) + if (.not. allocated(uinps)) allocate (uinps(3,maxatm)) c c find and store the atomic dipole polarizability parameters c @@ -245,14 +245,14 @@ subroutine polargrp maxp12 = 50 maxp13 = 50 maxp14 = 50 - if (.not. allocated(np11)) allocate (np11(n)) - if (.not. allocated(np12)) allocate (np12(n)) - if (.not. allocated(np13)) allocate (np13(n)) - if (.not. allocated(np14)) allocate (np14(n)) - if (.not. allocated(ip11)) allocate (ip11(maxp11,n)) - if (.not. allocated(ip12)) allocate (ip12(maxp12,n)) - if (.not. allocated(ip13)) allocate (ip13(maxp13,n)) - if (.not. allocated(ip14)) allocate (ip14(maxp14,n)) + if (.not. allocated(np11)) allocate (np11(maxatm)) + if (.not. allocated(np12)) allocate (np12(maxatm)) + if (.not. allocated(np13)) allocate (np13(maxatm)) + if (.not. allocated(np14)) allocate (np14(maxatm)) + if (.not. allocated(ip11)) allocate (ip11(maxp11,maxatm)) + if (.not. allocated(ip12)) allocate (ip12(maxp12,maxatm)) + if (.not. allocated(ip13)) allocate (ip13(maxp13,maxatm)) + if (.not. allocated(ip14)) allocate (ip14(maxp14,maxatm)) c c find the directly connected group members for each atom c diff --git a/source/ksolv.f b/source/ksolv.f index 55fe261fd..97c639c74 100644 --- a/source/ksolv.f +++ b/source/ksolv.f @@ -215,11 +215,6 @@ subroutine ksa integer atmnum c c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(rsolv)) allocate (rsolv(n)) - if (.not. allocated(asolv)) allocate (asolv(n)) -c c assign the Eisenberg-McLachlan ASP solvation parameters; c parameters only available for protein-peptide groups c @@ -373,23 +368,6 @@ subroutine kgb logical amide c c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(rsolv)) allocate (rsolv(n)) - if (.not. allocated(asolv)) allocate (asolv(n)) - if (.not. allocated(rborn)) allocate (rborn(n)) - if (.not. allocated(drb)) allocate (drb(n)) - if (.not. allocated(drobc)) allocate (drobc(n)) - if (.not. allocated(gpol)) allocate (gpol(n)) - if (.not. allocated(shct)) allocate (shct(n)) - if (.not. allocated(aobc)) allocate (aobc(n)) - if (.not. allocated(bobc)) allocate (bobc(n)) - if (.not. allocated(gobc)) allocate (gobc(n)) - if (.not. allocated(vsolv)) allocate (vsolv(n)) - if (.not. allocated(wace)) allocate (wace(maxclass,maxclass)) - if (.not. allocated(s2ace)) allocate (s2ace(maxclass,maxclass)) - if (.not. allocated(uace)) allocate (uace(maxclass,maxclass)) -c c set offset and scaling values for analytical Still method c if (borntyp .eq. 'STILL') then @@ -823,14 +801,6 @@ subroutine kgk character*120 string c c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(rsolv)) allocate (rsolv(n)) - if (.not. allocated(rborn)) allocate (rborn(n)) - if (.not. allocated(drb)) allocate (drb(n)) - if (.not. allocated(drbp)) allocate (drbp(n)) - if (.not. allocated(shct)) allocate (shct(n)) -c c set default value for exponent in the GB/GK function c gkc = 2.455d0 @@ -1134,8 +1104,8 @@ subroutine kgk c subroutine kpb use sizes - use atomid use atoms + use atomid use bath use couple use gkstuf @@ -1172,11 +1142,6 @@ subroutine kpb character*120 string c c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(rsolv)) allocate (rsolv(n)) - if (.not. allocated(shct)) allocate (shct(n)) -c c assign some default APBS configuration parameters c pbtyp = 'LPBE' @@ -1736,19 +1701,18 @@ subroutine knp stcut = cross + 3.9d0 stoff = cross - 3.5d0 c -c perform dynamic allocation of some global arrays -c - if (.not. allocated(asolv)) allocate (asolv(n)) - if (.not. allocated(rcav)) allocate (rcav(n)) - if (.not. allocated(rdisp)) allocate (rdisp(n)) - if (.not. allocated(cdisp)) allocate (cdisp(n)) -c c assign surface area factors for nonpolar solvation c do i = 1, n asolv(i) = surften end do c +c perform dynamic allocation of some global arrays +c + if (.not. allocated(rcav)) allocate (rcav(n)) + if (.not. allocated(rdisp)) allocate (rdisp(n)) + if (.not. allocated(cdisp)) allocate (cdisp(n)) +c c set cavity and dispersion radii for nonpolar solvation c do i = 1, n diff --git a/source/kstrbnd.f b/source/kstrbnd.f index 60c496c69..c063db5bd 100644 --- a/source/kstrbnd.f +++ b/source/kstrbnd.f @@ -112,8 +112,8 @@ subroutine kstrbnd c c perform dynamic allocation of some global arrays c - if (.not. allocated(isb)) allocate (isb(3,nangle)) - if (.not. allocated(sbk)) allocate (sbk(2,nangle)) + if (.not. allocated(isb)) allocate (isb(3,maxang)) + if (.not. allocated(sbk)) allocate (sbk(2,maxang)) c c use special stretch-bend parameter assignment method for MMFF c diff --git a/source/kstrtor.f b/source/kstrtor.f index e5af93b5c..31c07f033 100644 --- a/source/kstrtor.f +++ b/source/kstrtor.f @@ -148,8 +148,8 @@ subroutine kstrtor c c perform dynamic allocation of some global arrays c - if (.not. allocated(ist)) allocate (ist(4,ntors)) - if (.not. allocated(kst)) allocate (kst(9,ntors)) + if (.not. allocated(ist)) allocate (ist(4,maxtors)) + if (.not. allocated(kst)) allocate (kst(9,maxtors)) c c assign the stretch-torsion parameters for each torsion c diff --git a/source/ktors.f b/source/ktors.f index f57e80cee..7536f7d74 100644 --- a/source/ktors.f +++ b/source/ktors.f @@ -198,12 +198,12 @@ subroutine ktors c c perform dynamic allocation of some global arrays c - if (.not. allocated(tors1)) allocate (tors1(4,ntors)) - if (.not. allocated(tors2)) allocate (tors2(4,ntors)) - if (.not. allocated(tors3)) allocate (tors3(4,ntors)) - if (.not. allocated(tors4)) allocate (tors4(4,ntors)) - if (.not. allocated(tors5)) allocate (tors5(4,ntors)) - if (.not. allocated(tors6)) allocate (tors6(4,ntors)) + if (.not. allocated(tors1)) allocate (tors1(4,maxtors)) + if (.not. allocated(tors2)) allocate (tors2(4,maxtors)) + if (.not. allocated(tors3)) allocate (tors3(4,maxtors)) + if (.not. allocated(tors4)) allocate (tors4(4,maxtors)) + if (.not. allocated(tors5)) allocate (tors5(4,maxtors)) + if (.not. allocated(tors6)) allocate (tors6(4,maxtors)) c c use special torsional parameter assignment method for MMFF c diff --git a/source/ktortor.f b/source/ktortor.f index 5e23d729e..fb8da5d04 100644 --- a/source/ktortor.f +++ b/source/ktortor.f @@ -143,7 +143,7 @@ subroutine ktortor c c perform dynamic allocation of some global arrays c - if (.not. allocated(itt)) allocate (itt(3,nbitor)) + if (.not. allocated(itt)) allocate (itt(3,maxbitor)) c c check whether each torsion-torsion parameter is periodic; c assumes the "tbf" array is sorted with both indices in diff --git a/source/kurey.f b/source/kurey.f index 778d1d261..c99af3c7e 100644 --- a/source/kurey.f +++ b/source/kurey.f @@ -105,9 +105,9 @@ subroutine kurey c c perform dynamic allocation of some global arrays c - if (.not. allocated(iury)) allocate (iury(3,nangle)) - if (.not. allocated(uk)) allocate (uk(nangle)) - if (.not. allocated(ul)) allocate (ul(nangle)) + if (.not. allocated(iury)) allocate (iury(3,maxang)) + if (.not. allocated(uk)) allocate (uk(maxang)) + if (.not. allocated(ul)) allocate (ul(maxang)) c c assign the Urey-Bradley parameters for each angle c diff --git a/source/kvdw.f b/source/kvdw.f index 42b42bcde..cceab8801 100644 --- a/source/kvdw.f +++ b/source/kvdw.f @@ -266,10 +266,10 @@ subroutine kvdw c c perform dynamic allocation of some global arrays c - if (.not. allocated(ivdw)) allocate (ivdw(n)) - if (.not. allocated(jvdw)) allocate (jvdw(n)) - if (.not. allocated(ired)) allocate (ired(n)) - if (.not. allocated(kred)) allocate (kred(n)) + if (.not. allocated(ivdw)) allocate (ivdw(maxatm)) + if (.not. allocated(jvdw)) allocate (jvdw(maxatm)) + if (.not. allocated(ired)) allocate (ired(maxatm)) + if (.not. allocated(kred)) allocate (kred(maxatm)) if (.not. allocated(radmin)) & allocate (radmin(maxclass,maxclass)) if (.not. allocated(epsilon)) diff --git a/source/lbfgs.f b/source/lbfgs.f index 0aa197cae..2868f61ac 100644 --- a/source/lbfgs.f +++ b/source/lbfgs.f @@ -44,6 +44,7 @@ c c subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) + use sizes use inform use iounit use keys @@ -88,6 +89,12 @@ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) c c initialize some values to be used below c + if (nvar .gt. maxvar) then + write (iout,10) + 10 format (/,' LBFGS -- Too many Parameters,', + & ' Increase the Value of MAXVAR') + return + end if ncalls = 0 rms = sqrt(dble(nvar)) if (coordtype .eq. 'CARTESIAN') then @@ -135,35 +142,35 @@ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) call upcase (keyword) string = record(next:120) if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then - read (string,*,err=10,end=10) msav + read (string,*,err=20,end=20) msav else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then msav = 0 else if (keyword(1:7) .eq. 'FCTMIN ') then - read (string,*,err=10,end=10) fctmin + read (string,*,err=20,end=20) fctmin else if (keyword(1:8) .eq. 'MAXITER ') then - read (string,*,err=10,end=10) maxiter + read (string,*,err=20,end=20) maxiter else if (keyword(1:8) .eq. 'STEPMAX ') then - read (string,*,err=10,end=10) stpmax + read (string,*,err=20,end=20) stpmax else if (keyword(1:8) .eq. 'STEPMIN ') then - read (string,*,err=10,end=10) stpmin + read (string,*,err=20,end=20) stpmin else if (keyword(1:6) .eq. 'CAPPA ') then - read (string,*,err=10,end=10) cappa + read (string,*,err=20,end=20) cappa else if (keyword(1:9) .eq. 'SLOPEMAX ') then - read (string,*,err=10,end=10) slpmax + read (string,*,err=20,end=20) slpmax else if (keyword(1:7) .eq. 'ANGMAX ') then - read (string,*,err=10,end=10) angmax + read (string,*,err=20,end=20) angmax else if (keyword(1:7) .eq. 'INTMAX ') then - read (string,*,err=10,end=10) intmax + read (string,*,err=20,end=20) intmax end if - 10 continue + 20 continue end do c c check the number of saved correction vectors c if (msav.lt.0 .or. msav.gt.min(nvar,maxsav)) then msav = min(nvar,maxsav) - write (iout,20) msav - 20 format (/,' LBFGS -- Number of Saved Vectors Set to', + write (iout,30) msav + 30 format (/,' LBFGS -- Number of Saved Vectors Set to', & ' the Maximum of',i5) end if c @@ -171,17 +178,17 @@ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) c if (iprint .gt. 0) then if (msav .eq. 0) then - write (iout,30) - 30 format (/,' Steepest Descent Gradient Optimization :') write (iout,40) - 40 format (/,' SD Iter F Value G RMS F Move', + 40 format (/,' Steepest Descent Gradient Optimization :') + write (iout,50) + 50 format (/,' SD Iter F Value G RMS F Move', & ' X Move Angle FG Call Comment',/) else - write (iout,50) - 50 format (/,' Limited Memory BFGS Quasi-Newton', - & ' Optimization :') write (iout,60) - 60 format (/,' QN Iter F Value G RMS F Move', + 60 format (/,' Limited Memory BFGS Quasi-Newton', + & ' Optimization :') + write (iout,70) + 70 format (/,' QN Iter F Value G RMS F Move', & ' X Move Angle FG Call Comment',/) end if end if @@ -221,11 +228,11 @@ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) c if (iprint .gt. 0) then if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then - write (iout,70) niter,f,g_rms,ncalls - 70 format (i6,f14.4,f11.4,29x,i7) - else write (iout,80) niter,f,g_rms,ncalls - 80 format (i6,d14.4,d11.4,29x,i7) + 80 format (i6,f14.4,f11.4,29x,i7) + else + write (iout,90) niter,f,g_rms,ncalls + 90 format (i6,d14.4,d11.4,29x,i7) end if end if c @@ -371,13 +378,13 @@ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) if (done .or. mod(niter,iprint).eq.0) then if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. & g_rms.lt.1.0d5 .and. f_move.lt.1.0d5) then - write (iout,90) niter,f,g_rms,f_move,x_move, + write (iout,100) niter,f,g_rms,f_move,x_move, & angle,ncalls,status - 90 format (i6,f14.4,f11.4,f11.4,f10.4,f8.2,i7,3x,a9) + 100 format (i6,f14.4,f11.4,f11.4,f10.4,f8.2,i7,3x,a9) else - write (iout,100) niter,f,g_rms,f_move,x_move, + write (iout,110) niter,f,g_rms,f_move,x_move, & angle,ncalls,status - 100 format (i6,d14.4,d11.4,d11.4,f10.4,f8.2,i7,3x,a9) + 110 format (i6,d14.4,d11.4,d11.4,f10.4,f8.2,i7,3x,a9) end if end if end if @@ -408,11 +415,11 @@ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave) minimum = f if (iprint .gt. 0) then if (status.eq.'SmallGrad' .or. status.eq.'SmallFct ') then - write (iout,110) status - 110 format (/,' LBFGS -- Normal Termination due to ',a9) - else write (iout,120) status - 120 format (/,' LBFGS -- Incomplete Convergence due to ',a9) + 120 format (/,' LBFGS -- Normal Termination due to ',a9) + else + write (iout,130) status + 130 format (/,' LBFGS -- Incomplete Convergence due to ',a9) end if end if return diff --git a/source/lights.f b/source/lights.f index 30e6526c8..75d179199 100644 --- a/source/lights.f +++ b/source/lights.f @@ -15,9 +15,6 @@ c "lights" computes the set of nearest neighbor interactions c using the method of lights algorithm c -c note this version generates each pair only once via setting -c of the negative x-coordinate boundaries -c c literature reference: c c F. Sullivan, R. D. Mountain and J. O'Connell, "Molecular @@ -33,9 +30,7 @@ subroutine lights (cutoff,nsite,xsort,ysort,zsort) use iounit use light implicit none - integer i,j,k - integer nsite - integer extent + integer i,j,k,nsite real*8 cutoff,box real*8 xcut,ycut,zcut real*8 xmove,ymove,zmove @@ -45,6 +40,9 @@ subroutine lights (cutoff,nsite,xsort,ysort,zsort) real*8, allocatable :: xfrac(:) real*8, allocatable :: yfrac(:) real*8, allocatable :: zfrac(:) + logical first + save first + data first / .true. / c c c check that maximum number of replicates is not exceeded @@ -174,33 +172,20 @@ subroutine lights (cutoff,nsite,xsort,ysort,zsort) c perform dynamic allocation of some global arrays c nlight = (ncell+1) * nsite - extent = 0 - if (allocated(rgx)) extent = size(rgx) - if (extent .lt. nlight) then - if (allocated(kbx)) deallocate (kbx) - if (allocated(kby)) deallocate (kby) - if (allocated(kbz)) deallocate (kbz) - if (allocated(kex)) deallocate (kex) - if (allocated(key)) deallocate (key) - if (allocated(kez)) deallocate (kez) - if (allocated(locx)) deallocate (locx) - if (allocated(locy)) deallocate (locy) - if (allocated(locz)) deallocate (locz) - if (allocated(rgx)) deallocate (rgx) - if (allocated(rgy)) deallocate (rgy) - if (allocated(rgz)) deallocate (rgz) - allocate (kbx(nsite)) - allocate (kby(nsite)) - allocate (kbz(nsite)) - allocate (kex(nsite)) - allocate (key(nsite)) - allocate (kez(nsite)) - allocate (locx(nlight)) - allocate (locy(nlight)) - allocate (locz(nlight)) - allocate (rgx(nlight)) - allocate (rgy(nlight)) - allocate (rgz(nlight)) + if (first) then + first = .false. + if (.not. allocated(kbx)) allocate (kbx(nsite)) + if (.not. allocated(kby)) allocate (kby(nsite)) + if (.not. allocated(kbz)) allocate (kbz(nsite)) + if (.not. allocated(kex)) allocate (kex(nsite)) + if (.not. allocated(key)) allocate (key(nsite)) + if (.not. allocated(kez)) allocate (kez(nsite)) + if (.not. allocated(locx)) allocate (locx(nlight)) + if (.not. allocated(locy)) allocate (locy(nlight)) + if (.not. allocated(locz)) allocate (locz(nlight)) + if (.not. allocated(rgx)) allocate (rgx(nlight)) + if (.not. allocated(rgy)) allocate (rgy(nlight)) + if (.not. allocated(rgz)) allocate (rgz(nlight)) end if c c sort the coordinate components into ascending order diff --git a/source/mdinit.f b/source/mdinit.f index bfe2ade02..898012ebf 100644 --- a/source/mdinit.f +++ b/source/mdinit.f @@ -40,12 +40,13 @@ subroutine mdinit use usage implicit none integer i,j,k,idyn - integer size,next + integer size,next,len_isok integer lext,freeunit real*8 e,ekt,qterm real*8 maxwell,speed real*8 vec(3) real*8, allocatable :: derivs(:,:) + real*8 tempstor,LkT,kT,random logical exist character*7 ext character*20 keyword @@ -68,6 +69,8 @@ subroutine mdinit use_pred = .false. polpred = 'LSQR' iprint = 100 + nrespa = 1 + len_nhc = 4 c c set default values for temperature and pressure control c @@ -147,6 +150,10 @@ subroutine mdinit call upcase (volscale) else if (keyword(1:9) .eq. 'PRINTOUT ') then read (string,*,err=10,end=10) iprint + else if (keyword(1:11) .eq. 'NHC-LENGTH ') then + read (string,*,err=10,end=10) len_nhc + else if (keyword(1:6) .eq. 'NRESPA ') then + read (string,*,err=10,end=10) nrespa end if 10 continue end do @@ -325,24 +332,6 @@ subroutine mdinit if (integrate .eq. 'GHMC') dorest = .false. if (isothermal .and. thermostat.eq.'ANDERSEN') dorest = .false. c -c perform dynamic allocation of some global arrays -c - if (integrate .eq. 'RIGIDBODY') then - if (.not. allocated(xcmo)) allocate (xcmo(n)) - if (.not. allocated(ycmo)) allocate (ycmo(n)) - if (.not. allocated(zcmo)) allocate (zcmo(n)) - if (.not. allocated(vcm)) allocate (vcm(3,ngrp)) - if (.not. allocated(wcm)) allocate (wcm(3,ngrp)) - if (.not. allocated(lm)) allocate (lm(3,ngrp)) - if (.not. allocated(vc)) allocate (vc(3,ngrp)) - if (.not. allocated(wc)) allocate (wc(3,ngrp)) - if (.not. allocated(linear)) allocate (linear(ngrp)) - else - if (.not. allocated(v)) allocate (v(3,n)) - if (.not. allocated(a)) allocate (a(3,n)) - if (.not. allocated(aalt)) allocate (aalt(3,n)) - end if -c c try to restart using prior velocities and accelerations c dynfile = filename(1:leng)//'.dyn' @@ -403,6 +392,110 @@ subroutine mdinit end do deallocate (derivs) if (nuse .eq. n) call mdrest (0) + +c +c set physical and extended velocities, +c physical accelerations, and extended +c masses for isokinetic dynamics +c + else if (integrate .eq. 'ISOKINETIC') then + kT = boltzmann*kelvin + len_isok=1 + LkT = len_isok*kT + allocate (derivs(3,n)) + call gradient (e,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + do k = 1, len_nhc + q_iso1(k,j,i) = kT*tautemp*tautemp + q_iso2(k,j,i) = kT*tautemp*tautemp + speed = maxwell (q_iso1(k,j,i),kelvin) + call ranvec (vec) + v_iso1(k,j,i) = speed * vec(1) + v_iso2(k,j,i) = speed * vec(2) + end do + + v(j,i) = random () + do k = 1, len_isok + v_iso1(k,j,i) = random () + end do + tempstor = v(j,i) * v(j,i) + do k = 1, len_isok + tempstor = tempstor + (v_iso1(k,j,i)*v_iso1(k,j,i)) + end do + tempstor = sqrt(tempstor) + v(j,i) = v(j,i) / tempstor + do k = 1, len_isok + v_iso1(k,j,i) = v_iso1(k,j,i) / tempstor + end do + v(j,i) = v(j,i) / (sqrt(mass(i)/LkT)) + do k = 1, len_isok + tempstor=(sqrt(q_iso1(k,j,i) + & /(kT*dble(len_isok+1)))) + v_iso1(k,j,i) = v_iso1(k,j,i)/tempstor + end do + tempstor = mass(i)*v(j,i)*v(j,i) + do k = 1, len_isok + tempstor = tempstor + dble(len_isok)*q_iso1(k,j,i)* + & v_iso1(k,j,i)*v_iso1(k,j,i)/dble(len_isok+1) + end do + do k = 1, len_nhc + end do + + end do + end do + deallocate (derivs) + call prtdyn + if (nuse .eq. n) call mdrest (0) + +c +c set physical and extended velocities, +c physical accelerations, and extended +c masses for stochastic isokinetic dynamics +c + else if (integrate .eq. 'STOCH-ISOK') then + kT = boltzmann*kelvin + LkT = len_nhc*kT + allocate (derivs(3,n)) + call gradient (e,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + do k = 1, len_nhc + q_iso1(k,j,i) = kT*tautemp*tautemp + q_iso2(k,j,i) = kT*tautemp*tautemp + v_iso2(k,j,i) = kT*random ()/q_iso2(k,j,i) + end do + v(j,i) = random () + do k = 1, len_nhc + v_iso1(k,j,i) = random () + end do + tempstor = v(j,i) * v(j,i) + do k = 1, len_nhc + tempstor = tempstor + (v_iso1(k,j,i)*v_iso1(k,j,i)) + end do + tempstor = sqrt(tempstor) + v(j,i) = v(j,i) / tempstor + do k = 1, len_nhc + v_iso1(k,j,i) = v_iso1(k,j,i) / tempstor + end do + v(j,i) = v(j,i) / (sqrt(mass(i)/LkT)) + do k = 1, len_nhc + tempstor=(sqrt(q_iso1(k,j,i)/(kT*dble(len_nhc+1)))) + v_iso1(k,j,i) = v_iso1(k,j,i)/tempstor + end do + tempstor = mass(i)*v(j,i)*v(j,i) + do k = 1, len_nhc + tempstor = tempstor + dble(len_nhc)*q_iso1(k,j,i)* + & v_iso1(k,j,i)*v_iso1(k,j,i)/dble(len_nhc+1) + end do + end do + end do + deallocate (derivs) + call prtdyn + if (nuse .eq. n) call mdrest (0) + c c set velocities and accelerations for Cartesian dynamics c diff --git a/source/mdstuf.f b/source/mdstuf.f index 909473907..1d8c8f634 100644 --- a/source/mdstuf.f +++ b/source/mdstuf.f @@ -20,13 +20,14 @@ c frcsave logical flag to save force vector components c uindsave logical flag to save induced atomic dipoles c integrate type of molecular dynamics integration algorithm -c +c nrespa number of small time steps per large time step in RESPA c module mdstuf implicit none integer nfree integer irest integer bmnmix + integer nrespa logical dorest logical velsave logical frcsave diff --git a/source/moldyn.f b/source/moldyn.f index c9b32fe61..bf0b2f74f 100644 --- a/source/moldyn.f +++ b/source/moldyn.f @@ -18,9 +18,10 @@ c c module moldyn + use sizes implicit none - real*8, allocatable :: v(:,:) - real*8, allocatable :: a(:,:) - real*8, allocatable :: aalt(:,:) + real*8 v(3,maxatm) + real*8 a(3,maxatm) + real*8 aalt(3,maxatm) save end diff --git a/source/nblist.f b/source/nblist.f index 8f3fee330..4e4140044 100644 --- a/source/nblist.f +++ b/source/nblist.f @@ -13,7 +13,7 @@ c c c "nblist" constructs and maintains nonbonded pair neighbor lists -c for vdw, electrostatic and polarization interactions +c for vdw and electrostatic interactions c c subroutine nblist @@ -51,6 +51,7 @@ subroutine vlist use boxes use iounit use neigh + use openmp use vdw implicit none integer i,j,k @@ -62,6 +63,7 @@ subroutine vlist real*8, allocatable :: xred(:) real*8, allocatable :: yred(:) real*8, allocatable :: zred(:) + logical parallel logical, allocatable :: update(:) c c @@ -94,11 +96,16 @@ subroutine vlist call fatal end if c +c choose between serial and parallel list building +c + parallel = .false. + if (nthread .gt. 1) parallel = .true. +c c perform a complete list build instead of an update c if (dovlst) then dovlst = .false. - if (octahedron) then + if (parallel .or. octahedron) then call vbuild (xred,yred,zred) else call vlight (xred,yred,zred) @@ -229,7 +236,7 @@ subroutine vbuild (xred,yred,zred) use neigh use vdw implicit none - integer i,k + integer i,j,k real*8 xi,yi,zi real*8 xr,yr,zr,r2 real*8 xred(*) @@ -239,7 +246,7 @@ subroutine vbuild (xred,yred,zred) c c store coordinates to reflect update of the site c -!$OMP PARALLEL default(shared) private(i,k,xi,yi,zi,xr,yr,zr,r2) +!$OMP PARALLEL default(shared) private(i,j,k,xi,yi,zi,xr,yr,zr,r2) !$OMP DO schedule(guided) do i = 1, nvdw xi = xred(i) @@ -251,7 +258,7 @@ subroutine vbuild (xred,yred,zred) c c generate all neighbors for the site being rebuilt c - nvlst(i) = 0 + j = 0 do k = i+1, nvdw xr = xi - xred(k) yr = yi - yred(k) @@ -259,10 +266,11 @@ subroutine vbuild (xred,yred,zred) call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. vbuf2) then - nvlst(i) = nvlst(i) + 1 - vlst(nvlst(i),i) = k + j = j + 1 + vlst(j,i) = k end if end do + nvlst(i) = j c c check to see if the neighbor list is too long c @@ -317,9 +325,10 @@ subroutine vlight (xred,yred,zred) c c perform dynamic allocation of some local arrays c - allocate (xsort(nvdw)) - allocate (ysort(nvdw)) - allocate (zsort(nvdw)) + nlight = (ncell+1) * nvdw + allocate (xsort(nlight)) + allocate (ysort(nlight)) + allocate (zsort(nlight)) c c transfer interaction site coordinates to sorting arrays c @@ -336,7 +345,7 @@ subroutine vlight (xred,yred,zred) c use the method of lights to generate neighbors c off = sqrt(vbuf2) - call lightn (off,nvdw,xsort,ysort,zsort) + call lights (off,nvdw,xsort,ysort,zsort) c c perform deallocation of some local arrays c @@ -346,9 +355,6 @@ subroutine vlight (xred,yred,zred) c c loop over all atoms computing the neighbor lists c -!$OMP PARALLEL default(shared) private(i,j,k,xi,yi,zi, -!$OMP& xr,yr,zr,r2,kgy,kgz,start,stop,repeat) -!$OMP DO schedule(guided) do i = 1, nvdw xi = xred(i) yi = yred(i) @@ -365,7 +371,6 @@ subroutine vlight (xred,yred,zred) 10 continue do j = start, stop k = locx(j) - if (k .le. i) goto 20 kgy = rgy(k) if (kby(i) .le. key(i)) then if (kgy.lt.kby(i) .or. kgy.gt.key(i)) goto 20 @@ -384,29 +389,34 @@ subroutine vlight (xred,yred,zred) call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. vbuf2) then - nvlst(i) = nvlst(i) + 1 - vlst(nvlst(i),i) = k + if (i .lt. k) then + nvlst(i) = nvlst(i) + 1 + vlst(nvlst(i),i) = k + else + nvlst(k) = nvlst(k) + 1 + vlst(nvlst(k),k) = i + end if end if 20 continue end do if (repeat) then repeat = .false. start = kbx(i) + 1 - stop = nvdw + stop = nlight goto 10 end if + end do c -c check to see if the neighbor list is too long +c check to see if the neighbor lists are too long c + do i = 1, nvdw if (nvlst(i) .ge. maxvlst) then write (iout,30) - 30 format (/,' VLIGHT -- Too many Neighbors;', + 30 format (/,' VFULL -- Too many Neighbors;', & ' Increase MAXVLST') call fatal end if end do -!$OMP END DO -!$OMP END PARALLEL return end c @@ -430,12 +440,14 @@ subroutine clist use charge use iounit use neigh + use openmp implicit none integer i,j,k integer ii,kk real*8 xi,yi,zi real*8 xr,yr,zr real*8 radius,r2 + logical parallel logical, allocatable :: update(:) c c @@ -454,11 +466,16 @@ subroutine clist call fatal end if c +c choose between serial and parallel list building +c + parallel = .false. + if (nthread .gt. 1) parallel = .true. +c c perform a complete list build instead of an update c if (doclst) then doclst = .false. - if (octahedron) then + if (parallel .or. octahedron) then call cbuild else call clight @@ -590,7 +607,7 @@ subroutine cbuild use iounit use neigh implicit none - integer i,k + integer i,j,k integer ii,kk real*8 xi,yi,zi real*8 xr,yr,zr,r2 @@ -598,7 +615,7 @@ subroutine cbuild c c store new coordinates to reflect update of the site c -!$OMP PARALLEL default(shared) private(i,k,ii,kk,xi,yi,zi,xr,yr,zr,r2) +!$OMP PARALLEL default(shared) private(i,j,k,ii,kk,xi,yi,zi,xr,yr,zr,r2) !$OMP DO schedule(guided) do i = 1, nion ii = kion(i) @@ -611,7 +628,7 @@ subroutine cbuild c c generate all neighbors for the site being rebuilt c - nelst(i) = 0 + j = 0 do k = i+1, nion kk = kion(k) xr = xi - x(kk) @@ -620,10 +637,11 @@ subroutine cbuild call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. cbuf2) then - nelst(i) = nelst(i) + 1 - elst(nelst(i),i) = k + j = j + 1 + elst(j,i) = k end if end do + nelst(i) = j c c check to see if the neighbor list is too long c @@ -676,9 +694,10 @@ subroutine clight c c perform dynamic allocation of some local arrays c - allocate (xsort(nion)) - allocate (ysort(nion)) - allocate (zsort(nion)) + nlight = (ncell+1) * nion + allocate (xsort(nlight)) + allocate (ysort(nlight)) + allocate (zsort(nlight)) c c transfer interaction site coordinates to sorting arrays c @@ -696,7 +715,7 @@ subroutine clight c use the method of lights to generate neighbors c off = sqrt(cbuf2) - call lightn (off,nion,xsort,ysort,zsort) + call lights (off,nion,xsort,ysort,zsort) c c perform deallocation of some local arrays c @@ -706,9 +725,6 @@ subroutine clight c c loop over all atoms computing the neighbor lists c -!$OMP PARALLEL default(shared) private(i,j,k,ii,kk,xi,yi,zi, -!$OMP& xr,yr,zr,r2,kgy,kgz,start,stop,repeat) -!$OMP DO schedule(guided) do i = 1, nion ii = kion(i) xi = x(ii) @@ -726,7 +742,6 @@ subroutine clight 10 continue do j = start, stop k = locx(j) - if (k .le. i) goto 20 kk = kion(k) kgy = rgy(k) if (kby(i) .le. key(i)) then @@ -746,29 +761,34 @@ subroutine clight call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. cbuf2) then - nelst(i) = nelst(i) + 1 - elst(nelst(i),i) = k + if (i .lt. k) then + nelst(i) = nelst(i) + 1 + elst(nelst(i),i) = k + else + nelst(k) = nelst(k) + 1 + elst(nelst(k),k) = i + end if end if 20 continue end do if (repeat) then repeat = .false. start = kbx(i) + 1 - stop = nion + stop = nlight goto 10 end if + end do c -c check to see if the neighbor list is too long +c check to see if the neighbor lists are too long c + do i = 1, nion if (nelst(i) .ge. maxelst) then write (iout,30) - 30 format (/,' CLIGHT -- Too many Neighbors;', + 30 format (/,' CFULL -- Too many Neighbors;', & ' Increase MAXELST') call fatal end if end do -!$OMP END DO -!$OMP END PARALLEL return end c @@ -792,12 +812,14 @@ subroutine mlist use iounit use mpole use neigh + use openmp implicit none integer i,j,k integer ii,kk real*8 xi,yi,zi real*8 xr,yr,zr real*8 radius,r2 + logical parallel logical, allocatable :: update(:) c c @@ -816,11 +838,16 @@ subroutine mlist call fatal end if c +c choose between serial and parallel list building +c + parallel = .false. + if (nthread .gt. 1) parallel = .true. +c c perform a complete list build instead of an update c if (domlst) then domlst = .false. - if (octahedron) then + if (parallel .or. octahedron) then call mbuild else call mlight @@ -952,15 +979,14 @@ subroutine mbuild use mpole use neigh implicit none - integer i,k - integer ii,kk + integer i,j,k,ii,kk real*8 xi,yi,zi real*8 xr,yr,zr,r2 c c c store new coordinates to reflect update of the site c -!$OMP PARALLEL default(shared) private(i,k,ii,kk,xi,yi,zi,xr,yr,zr,r2) +!$OMP PARALLEL default(shared) private(i,j,k,ii,kk,xi,yi,zi,xr,yr,zr,r2) !$OMP DO schedule(guided) do i = 1, npole ii = ipole(i) @@ -973,7 +999,7 @@ subroutine mbuild c c generate all neighbors for the site being rebuilt c - nelst(i) = 0 + j = 0 do k = i+1, npole kk = ipole(k) xr = xi - x(kk) @@ -982,10 +1008,11 @@ subroutine mbuild call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. mbuf2) then - nelst(i) = nelst(i) + 1 - elst(nelst(i),i) = k + j = j + 1 + elst(j,i) = k end if end do + nelst(i) = j c c check to see if the neighbor list is too long c @@ -1038,9 +1065,10 @@ subroutine mlight c c perform dynamic allocation of some local arrays c - allocate (xsort(npole)) - allocate (ysort(npole)) - allocate (zsort(npole)) + nlight = (ncell+1) * npole + allocate (xsort(nlight)) + allocate (ysort(nlight)) + allocate (zsort(nlight)) c c transfer interaction site coordinates to sorting arrays c @@ -1058,7 +1086,7 @@ subroutine mlight c use the method of lights to generate neighbors c off = sqrt(mbuf2) - call lightn (off,npole,xsort,ysort,zsort) + call lights (off,npole,xsort,ysort,zsort) c c perform deallocation of some local arrays c @@ -1068,9 +1096,6 @@ subroutine mlight c c loop over all atoms computing the neighbor lists c -!$OMP PARALLEL default(shared) private(i,j,k,ii,kk,xi,yi,zi, -!$OMP& xr,yr,zr,r2,kgy,kgz,start,stop,repeat) -!$OMP DO schedule(guided) do i = 1, npole ii = ipole(i) xi = x(ii) @@ -1088,7 +1113,6 @@ subroutine mlight 10 continue do j = start, stop k = locx(j) - if (k .le. i) goto 20 kk = ipole(k) kgy = rgy(k) if (kby(i) .le. key(i)) then @@ -1108,29 +1132,34 @@ subroutine mlight call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. mbuf2) then - nelst(i) = nelst(i) + 1 - elst(nelst(i),i) = k + if (i .lt. k) then + nelst(i) = nelst(i) + 1 + elst(nelst(i),i) = k + else + nelst(k) = nelst(k) + 1 + elst(nelst(k),k) = i + end if end if 20 continue end do if (repeat) then repeat = .false. start = kbx(i) + 1 - stop = npole + stop = nlight goto 10 end if + end do c -c check to see if the neighbor list is too long +c check to see if the neighbor lists are too long c + do i = 1, npole if (nelst(i) .ge. maxelst) then write (iout,30) - 30 format (/,' MLIGHT -- Too many Neighbors;', + 30 format (/,' MFULL -- Too many Neighbors;', & ' Increase MAXELST') call fatal end if end do -!$OMP END DO -!$OMP END PARALLEL return end c @@ -1154,12 +1183,14 @@ subroutine ulist use iounit use mpole use neigh + use openmp implicit none integer i,j,k integer ii,kk real*8 xi,yi,zi real*8 xr,yr,zr real*8 radius,r2 + logical parallel logical, allocatable :: update(:) c c @@ -1178,11 +1209,16 @@ subroutine ulist call fatal end if c +c choose between serial and parallel list building +c + parallel = .false. + if (nthread .gt. 1) parallel = .true. +c c perform a complete list build instead of an update c if (doulst) then doulst = .false. - if (octahedron) then + if (parallel .or. octahedron) then call ubuild else call ulight @@ -1314,15 +1350,14 @@ subroutine ubuild use mpole use neigh implicit none - integer i,k - integer ii,kk + integer i,j,k,ii,kk real*8 xi,yi,zi real*8 xr,yr,zr,r2 c c c store new coordinates to reflect update of the site c -!$OMP PARALLEL default(shared) private(i,k,ii,kk,xi,yi,zi,xr,yr,zr,r2) +!$OMP PARALLEL default(shared) private(i,j,k,ii,kk,xi,yi,zi,xr,yr,zr,r2) !$OMP DO schedule(guided) do i = 1, npole ii = ipole(i) @@ -1335,7 +1370,7 @@ subroutine ubuild c c generate all neighbors for the site being rebuilt c - nulst(i) = 0 + j = 0 do k = i+1, npole kk = ipole(k) xr = xi - x(kk) @@ -1344,10 +1379,11 @@ subroutine ubuild call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. ubuf2) then - nulst(i) = nulst(i) + 1 - ulst(nulst(i),i) = k + j = j + 1 + ulst(j,i) = k end if end do + nulst(i) = j c c check to see if the neighbor list is too long c @@ -1401,9 +1437,10 @@ subroutine ulight c c perform dynamic allocation of some local arrays c - allocate (xsort(npole)) - allocate (ysort(npole)) - allocate (zsort(npole)) + nlight = (ncell+1) * npole + allocate (xsort(nlight)) + allocate (ysort(nlight)) + allocate (zsort(nlight)) c c transfer interaction site coordinates to sorting arrays c @@ -1421,7 +1458,7 @@ subroutine ulight c use the method of lights to generate neighbors c off = sqrt(ubuf2) - call lightn (off,npole,xsort,ysort,zsort) + call lights (off,npole,xsort,ysort,zsort) c c perform deallocation of some local arrays c @@ -1431,9 +1468,6 @@ subroutine ulight c c loop over all atoms computing the neighbor lists c -!$OMP PARALLEL default(shared) private(i,j,k,ii,kk,xi,yi,zi, -!$OMP& xr,yr,zr,r2,kgy,kgz,start,stop,repeat) -!$OMP DO schedule(guided) do i = 1, npole ii = ipole(i) xi = x(ii) @@ -1451,7 +1485,6 @@ subroutine ulight 10 continue do j = start, stop k = locx(j) - if (k .le. i) goto 20 kk = ipole(k) kgy = rgy(k) if (kby(i) .le. key(i)) then @@ -1471,355 +1504,34 @@ subroutine ulight call imagen (xr,yr,zr) r2 = xr*xr + yr*yr + zr*zr if (r2 .le. ubuf2) then - nulst(i) = nulst(i) + 1 - ulst(nulst(i),i) = k + if (i .lt. k) then + nulst(i) = nulst(i) + 1 + ulst(nulst(i),i) = k + else + nulst(k) = nulst(k) + 1 + ulst(nulst(k),k) = i + end if end if 20 continue end do if (repeat) then repeat = .false. start = kbx(i) + 1 - stop = npole + stop = nlight goto 10 end if + end do c -c check to see if the neighbor list is too long +c check to see if the neighbor lists are too long c + do i = 1, npole if (nulst(i) .ge. maxulst) then write (iout,30) - 30 format (/,' ULIGHT -- Too many Neighbors;', + 30 format (/,' UFULL -- Too many Neighbors;', & ' Increase MAXULST') call fatal end if end do -!$OMP END DO -!$OMP END PARALLEL - return - end -c -c -c ################################################################# -c ## ## -c ## subroutine lightn -- method of lights for list building ## -c ## ## -c ################################################################# -c -c -c "lightn" computes the set of nearest neighbor interactions -c using the method of lights algorithm -c -c note this is a special version for neighbor list generation -c which generates each pair in both directions, (A,B) and (B,A) -c -c literature reference: -c -c F. Sullivan, R. D. Mountain and J. O'Connell, "Molecular -c Dynamics on Vector Computers", Journal of Computational -c Physics, 61, 138-153 (1985) -c -c - subroutine lightn (cutoff,nsite,xsort,ysort,zsort) - use sizes - use bound - use boxes - use cell - use iounit - use light - implicit none - integer i,j,k - integer nsite - integer extent - real*8 cutoff,box - real*8 xcut,ycut,zcut - real*8 xsort(*) - real*8 ysort(*) - real*8 zsort(*) - real*8, allocatable :: xfrac(:) - real*8, allocatable :: yfrac(:) - real*8, allocatable :: zfrac(:) -c -c -c truncated octahedron periodicity is not handled at present -c - if (use_bounds) then - if (octahedron) then - write (iout,10) - 10 format (/,' LIGHTS -- Truncated Octahedron not', - & ' Supported by Method of Lights') - call fatal - end if - end if -c -c set the light width based on input distance cutoff -c - xcut = cutoff - ycut = cutoff - zcut = cutoff - if (use_bounds) then - if (monoclinic) then - zcut = zcut / beta_sin - xcut = xcut + zcut*abs(beta_cos) - else if (triclinic) then - zcut = zcut / gamma_term - ycut = (ycut + zcut*abs(beta_term)) / gamma_sin - xcut = xcut + ycut*abs(gamma_cos) + zcut*abs(beta_cos) - end if - xcut = min(xcut,xcell2) - ycut = min(ycut,ycell2) - zcut = min(zcut,zcell2) - end if -c -c perform dynamic allocation of some local arrays -c - allocate (xfrac(nsite)) - allocate (yfrac(nsite)) - allocate (zfrac(nsite)) -c -c find fractional coordinates for the unitcell atoms -c - if (use_bounds) then - if (orthogonal) then - do i = 1, nsite - zfrac(i) = zsort(i) - yfrac(i) = ysort(i) - xfrac(i) = xsort(i) - end do - else if (monoclinic) then - do i = 1, nsite - zfrac(i) = zsort(i) / beta_sin - yfrac(i) = ysort(i) - xfrac(i) = xsort(i) - zfrac(i)*beta_cos - end do - else if (triclinic) then - do i = 1, nsite - zfrac(i) = zsort(i) / gamma_term - yfrac(i) = (ysort(i) - zfrac(i)*beta_term) / gamma_sin - xfrac(i) = xsort(i) - yfrac(i)*gamma_cos - & - zfrac(i)*beta_cos - end do - end if - end if -c -c use images to move coordinates into periodic cell -c - if (use_bounds) then - do i = 1, nsite - xsort(i) = xfrac(i) - ysort(i) = yfrac(i) - zsort(i) = zfrac(i) - do while (abs(xsort(i)) .gt. xcell2) - xsort(i) = xsort(i) - sign(xcell,xsort(i)) - end do - do while (abs(ysort(i)) .gt. ycell2) - ysort(i) = ysort(i) - sign(ycell,ysort(i)) - end do - do while (abs(zsort(i)) .gt. zcell2) - zsort(i) = zsort(i) - sign(zcell,zsort(i)) - end do - end do - end if -c -c perform deallocation of some local arrays -c - deallocate (xfrac) - deallocate (yfrac) - deallocate (zfrac) -c -c perform dynamic allocation of some global arrays -c - nlight = nsite - extent = 0 - if (allocated(rgx)) extent = size(rgx) - if (extent .lt. nlight) then - if (allocated(kbx)) deallocate (kbx) - if (allocated(kby)) deallocate (kby) - if (allocated(kbz)) deallocate (kbz) - if (allocated(kex)) deallocate (kex) - if (allocated(key)) deallocate (key) - if (allocated(kez)) deallocate (kez) - if (allocated(locx)) deallocate (locx) - if (allocated(locy)) deallocate (locy) - if (allocated(locz)) deallocate (locz) - if (allocated(rgx)) deallocate (rgx) - if (allocated(rgy)) deallocate (rgy) - if (allocated(rgz)) deallocate (rgz) - allocate (kbx(nsite)) - allocate (kby(nsite)) - allocate (kbz(nsite)) - allocate (kex(nsite)) - allocate (key(nsite)) - allocate (kez(nsite)) - allocate (locx(nlight)) - allocate (locy(nlight)) - allocate (locz(nlight)) - allocate (rgx(nlight)) - allocate (rgy(nlight)) - allocate (rgz(nlight)) - end if -c -c sort the coordinate components into ascending order -c - call sort2 (nlight,xsort,locx) - call sort2 (nlight,ysort,locy) - call sort2 (nlight,zsort,locz) -c -c index the position of each atom in the sorted coordinates -c - do i = 1, nlight - rgx(locx(i)) = i - rgy(locy(i)) = i - rgz(locz(i)) = i - end do -c -c find the negative x-coordinate boundary for each atom -c - j = nlight - box = 0.0d0 - do i = nlight, 1, -1 - k = locx(i) - do while (xsort(i)-xsort(j)+box .le. xcut) - if (j .eq. 1) then - if (use_bounds) then - j = nlight + 1 - box = xcell - end if - end if - j = j - 1 - if (j .lt. 1) goto 20 - end do - 20 continue - j = j + 1 - if (j .gt. nlight) then - j = 1 - box = 0.0d0 - end if - kbx(k) = j - end do -c -c find the positive x-coordinate boundary for each atom -c - j = 1 - box = 0.0d0 - do i = 1, nlight - k = locx(i) - do while (xsort(j)-xsort(i)+box .lt. xcut) - if (j .eq. nlight) then - if (use_bounds) then - j = 0 - box = xcell - end if - end if - j = j + 1 - if (j .gt. nlight) goto 30 - end do - 30 continue - j = j - 1 - if (j .lt. 1) then - j = nlight - box = 0.0d0 - end if - kex(k) = j - end do -c -c find the negative y-coordinate boundary for each atom -c - j = nlight - box = 0.0d0 - do i = nlight, 1, -1 - k = locy(i) - do while (ysort(i)-ysort(j)+box .le. ycut) - if (j .eq. 1) then - if (use_bounds) then - j = nlight + 1 - box = ycell - end if - end if - j = j - 1 - if (j .lt. 1) goto 40 - end do - 40 continue - j = j + 1 - if (j .gt. nlight) then - j = 1 - box = 0.0d0 - end if - kby(k) = j - end do -c -c find the positive y-coordinate boundary for each atom -c - j = 1 - box = 0.0d0 - do i = 1, nlight - k = locy(i) - do while (ysort(j)-ysort(i)+box .lt. ycut) - if (j .eq. nlight) then - if (use_bounds) then - j = 0 - box = ycell - end if - end if - j = j + 1 - if (j .gt. nlight) goto 50 - end do - 50 continue - j = j - 1 - if (j .lt. 1) then - j = nlight - box = 0.0d0 - end if - key(k) = j - end do -c -c find the negative z-coordinate boundary for each atom -c - j = nlight - box = 0.0d0 - do i = nlight, 1, -1 - k = locz(i) - do while (zsort(i)-zsort(j)+box .le. zcut) - if (j .eq. 1) then - if (use_bounds) then - j = nlight + 1 - box = zcell - end if - end if - j = j - 1 - if (j .lt. 1) goto 60 - end do - 60 continue - j = j + 1 - if (j .gt. nlight) then - j = 1 - box = 0.0d0 - end if - kbz(k) = j - end do -c -c find the positive z-coordinate boundary for each atom -c - j = 1 - box = 0.0d0 - do i = 1, nlight - k = locz(i) - do while (zsort(j)-zsort(i)+box .lt. zcut) - if (j .eq. nlight) then - if (use_bounds) then - j = 0 - box = zcell - end if - end if - j = j + 1 - if (j .gt. nlight) goto 70 - end do - 70 continue - j = j - 1 - if (j .lt. 1) then - j = nlight - box = 0.0d0 - end if - kez(k) = j - end do return end c @@ -1833,9 +1545,7 @@ subroutine lightn (cutoff,nsite,xsort,ysort,zsort) c c "imagen" takes the components of pairwise distance between c two points and converts to the components of the minimum -c image distance -c -c note this is a fast version for neighbor list generation +c image distance; fast version for neighbor list generation c which only returns the correct component magnitudes c c diff --git a/source/nucleic.f b/source/nucleic.f index b23ad5fe8..a3c19c4c3 100644 --- a/source/nucleic.f +++ b/source/nucleic.f @@ -131,6 +131,7 @@ program nucleic c c subroutine getseqn + use sizes use iounit use nucleo use resdue diff --git a/source/params.f b/source/params.f index 5c9d48d2c..f666bb654 100644 --- a/source/params.f +++ b/source/params.f @@ -12,16 +12,13 @@ c ############################################################## c c -c maxprm maximum number of lines in the parameter file -c c nprm number of nonblank lines in the parameter file c prmline contents of each individual parameter file line c c module params + use sizes implicit none - integer maxprm - parameter (maxprm=25000) integer nprm character*120 prmline(maxprm) save diff --git a/source/protein.f b/source/protein.f index 2f82b95c7..70e3aa342 100644 --- a/source/protein.f +++ b/source/protein.f @@ -140,6 +140,7 @@ program protein c c subroutine getseq + use sizes use iounit use phipsi use resdue diff --git a/source/ptable.f b/source/ptable.f index 593b18d24..e127c793a 100644 --- a/source/ptable.f +++ b/source/ptable.f @@ -12,15 +12,12 @@ c ############################################################### c c -c maxele maximum number of elements from periodic table -c c elemnt atomic symbol for each chemical element c c module ptable + use sizes implicit none - integer maxele - parameter (maxele=112) character*3 elemnt(maxele) save end diff --git a/source/readgau.f b/source/readgau.f index 8e68fb1c0..d5672ff20 100644 --- a/source/readgau.f +++ b/source/readgau.f @@ -33,7 +33,6 @@ subroutine readgau logical hasinputxyz logical hasmp2 logical exist - real*8 xtmp,ytmp,ztmp real*8 frcunit,hessunit character*4 arcstart character*6 gname @@ -105,7 +104,7 @@ subroutine readgau read (igau,50,err=70,end=70) record 50 format (a120) read (record,*,err=60,end=60) itmp,jtmp,ktmp, - & xtmp,ytmp,ztmp + & gx(i),gy(i),gz(i) if (jtmp .le. 0) goto 60 i = i + 1 end do @@ -147,7 +146,7 @@ subroutine readgau read (igau,100,err=220,end=220) record 100 format (a120) read (record,*,err=110,end=110) itmp,jtmp,ktmp, - & gx(i),gy(i),gz(i) + & gx(i),gy(i),gz(i) if (jtmp .le. 0) goto 110 i = i + 1 end do diff --git a/source/readpdb.f b/source/readpdb.f index 6cc7006f8..24569d688 100644 --- a/source/readpdb.f +++ b/source/readpdb.f @@ -17,6 +17,7 @@ c c subroutine readpdb (ipdb) + use sizes use files use inform use iounit @@ -375,6 +376,7 @@ subroutine readpdb (ipdb) c c subroutine scanpdb (ipdb) + use sizes use iounit use pdb use sequen @@ -594,6 +596,7 @@ subroutine scanpdb (ipdb) c c subroutine fixpdb (resname,atmname) + use sizes use resdue implicit none integer i diff --git a/source/readseq.f b/source/readseq.f index 4694c0354..93d7ce169 100644 --- a/source/readseq.f +++ b/source/readseq.f @@ -19,6 +19,7 @@ c c subroutine readseq (iseq) + use sizes use files use iounit use resdue diff --git a/source/resdue.f b/source/resdue.f index 969a3dc24..c546035fc 100644 --- a/source/resdue.f +++ b/source/resdue.f @@ -12,60 +12,54 @@ c ################################################################ c c -c maxamino maximum number of amino acid residue types -c maxnuc maximum number of nucleic acid residue types -c -c ntyp biotypes for mid-chain peptide backbone N atoms -c catyp biotypes for mid-chain peptide backbone CA atoms -c ctyp biotypes for mid-chain peptide backbone C atoms -c hntyp biotypes for mid-chain peptide backbone HN atoms -c otyp biotypes for mid-chain peptide backbone O atoms -c hatyp biotypes for mid-chain peptide backbone HA atoms -c cbtyp biotypes for mid-chain peptide backbone CB atoms -c nntyp biotypes for N-terminal peptide backbone N atoms -c cantyp biotypes for N-terminal peptide backbone CA atoms -c cntyp biotypes for N-terminal peptide backbone C atoms -c hnntyp biotypes for N-terminal peptide backbone HN atoms -c ontyp biotypes for N-terminal peptide backbone O atoms -c hantyp biotypes for N-terminal peptide backbone HA atoms -c nctyp biotypes for C-terminal peptide backbone N atoms -c cactyp biotypes for C-terminal peptide backbone CA atoms -c cctyp biotypes for C-terminal peptide backbone C atoms -c hnctyp biotypes for C-terminal peptide backbone HN atoms -c octyp biotypes for C-terminal peptide backbone O atoms -c hactyp biotypes for C-terminal peptide backbone HA atoms -c o5typ biotypes for nucleotide backbone and sugar O5' atoms -c c5typ biotypes for nucleotide backbone and sugar C5' atoms -c h51typ biotypes for nucleotide backbone and sugar H5' atoms -c h52typ biotypes for nucleotide backbone and sugar H5'' atoms -c c4typ biotypes for nucleotide backbone and sugar C4' atoms -c h4typ biotypes for nucleotide backbone and sugar H4' atoms -c o4typ biotypes for nucleotide backbone and sugar O4' atoms -c c1typ biotypes for nucleotide backbone and sugar C1' atoms -c h1typ biotypes for nucleotide backbone and sugar H1' atoms -c c3typ biotypes for nucleotide backbone and sugar C3' atoms -c h3typ biotypes for nucleotide backbone and sugar H3' atoms -c c2typ biotypes for nucleotide backbone and sugar C2' atoms -c h21typ biotypes for nucleotide backbone and sugar H2' atoms -c o2typ biotypes for nucleotide backbone and sugar O2' atoms -c h22typ biotypes for nucleotide backbone and sugar H2'' atoms -c o3typ biotypes for nucleotide backbone and sugar O3' atoms -c ptyp biotypes for nucleotide backbone and sugar P atoms -c optyp biotypes for nucleotide backbone and sugar OP atoms -c h5ttyp biotypes for nucleotide backbone and sugar H5T atoms -c h3ttyp biotypes for nucleotide backbone and sugar H3T atoms -c amino three-letter abbreviations for amino acids types -c nuclz three-letter abbreviations for nucleic acids types -c amino1 one-letter abbreviations for amino acids types -c nuclz1 one-letter abbreviations for nucleic acids types +c ntyp biotypes for mid-chain peptide backbone N atoms +c catyp biotypes for mid-chain peptide backbone CA atoms +c ctyp biotypes for mid-chain peptide backbone C atoms +c hntyp biotypes for mid-chain peptide backbone HN atoms +c otyp biotypes for mid-chain peptide backbone O atoms +c hatyp biotypes for mid-chain peptide backbone HA atoms +c cbtyp biotypes for mid-chain peptide backbone CB atoms +c nntyp biotypes for N-terminal peptide backbone N atoms +c cantyp biotypes for N-terminal peptide backbone CA atoms +c cntyp biotypes for N-terminal peptide backbone C atoms +c hnntyp biotypes for N-terminal peptide backbone HN atoms +c ontyp biotypes for N-terminal peptide backbone O atoms +c hantyp biotypes for N-terminal peptide backbone HA atoms +c nctyp biotypes for C-terminal peptide backbone N atoms +c cactyp biotypes for C-terminal peptide backbone CA atoms +c cctyp biotypes for C-terminal peptide backbone C atoms +c hnctyp biotypes for C-terminal peptide backbone HN atoms +c octyp biotypes for C-terminal peptide backbone O atoms +c hactyp biotypes for C-terminal peptide backbone HA atoms +c o5typ biotypes for nucleotide backbone and sugar O5' atoms +c c5typ biotypes for nucleotide backbone and sugar C5' atoms +c h51typ biotypes for nucleotide backbone and sugar H5' atoms +c h52typ biotypes for nucleotide backbone and sugar H5'' atoms +c c4typ biotypes for nucleotide backbone and sugar C4' atoms +c h4typ biotypes for nucleotide backbone and sugar H4' atoms +c o4typ biotypes for nucleotide backbone and sugar O4' atoms +c c1typ biotypes for nucleotide backbone and sugar C1' atoms +c h1typ biotypes for nucleotide backbone and sugar H1' atoms +c c3typ biotypes for nucleotide backbone and sugar C3' atoms +c h3typ biotypes for nucleotide backbone and sugar H3' atoms +c c2typ biotypes for nucleotide backbone and sugar C2' atoms +c h21typ biotypes for nucleotide backbone and sugar H2' atoms +c o2typ biotypes for nucleotide backbone and sugar O2' atoms +c h22typ biotypes for nucleotide backbone and sugar H2'' atoms +c o3typ biotypes for nucleotide backbone and sugar O3' atoms +c ptyp biotypes for nucleotide backbone and sugar P atoms +c optyp biotypes for nucleotide backbone and sugar OP atoms +c h5ttyp biotypes for nucleotide backbone and sugar H5T atoms +c h3ttyp biotypes for nucleotide backbone and sugar H3T atoms +c amino three-letter abbreviations for amino acids types +c nuclz three-letter abbreviations for nucleic acids types +c amino1 one-letter abbreviations for amino acids types +c nuclz1 one-letter abbreviations for nucleic acids types c c module resdue + use sizes implicit none - integer maxamino - integer maxnuc - parameter (maxamino=38) - parameter (maxnuc=12) integer ntyp(maxamino) integer catyp(maxamino) integer ctyp(maxamino) diff --git a/source/rgddyn.f b/source/rgddyn.f index c24fcd7a1..75d738bff 100644 --- a/source/rgddyn.f +++ b/source/rgddyn.f @@ -24,15 +24,16 @@ c c module rgddyn + use sizes implicit none - real*8, allocatable :: xcmo(:) - real*8, allocatable :: ycmo(:) - real*8, allocatable :: zcmo(:) - real*8, allocatable :: vcm(:,:) - real*8, allocatable :: wcm(:,:) - real*8, allocatable :: lm(:,:) - real*8, allocatable :: vc(:,:) - real*8, allocatable :: wc(:,:) - logical, allocatable :: linear(:) + real*8 xcmo(maxatm) + real*8 ycmo(maxatm) + real*8 zcmo(maxatm) + real*8 vcm(3,maxgrp) + real*8 wcm(3,maxgrp) + real*8 lm(3,maxgrp) + real*8 vc(3,maxgrp) + real*8 wc(3,maxgrp) + logical linear(maxgrp) save end diff --git a/source/rings.f b/source/rings.f index e59654aa2..9880d9105 100644 --- a/source/rings.f +++ b/source/rings.f @@ -21,6 +21,7 @@ c c subroutine rings + use sizes use angbnd use atoms use bitor @@ -36,7 +37,6 @@ subroutine rings integer id,ie,ig integer list1,list2 integer list3,list4 - integer maxring integer, allocatable :: list(:) logical reduce c @@ -58,7 +58,6 @@ subroutine rings c c perform dynamic allocation of some global arrays c - maxring = 10000 if (.not. allocated(iring3)) allocate (iring3(3,maxring)) if (.not. allocated(iring4)) allocate (iring4(4,maxring)) if (.not. allocated(iring5)) allocate (iring5(5,maxring)) @@ -76,8 +75,7 @@ subroutine rings nring3 = nring3 + 1 if (nring3 .gt. maxring) then write (iout,10) - 10 format (/,' RINGS -- Too many 3-Membered Rings;', - & ' Increase MAXRING') + 10 format (/,' RINGS -- Too many 3-Membered Rings') call fatal end if iring3(1,nring3) = ia @@ -110,8 +108,7 @@ subroutine rings nring4 = nring4 + 1 if (nring4 .gt. maxring) then write (iout,30) - 30 format (/,' RINGS -- Too many 4-Membered Rings;', - & ' Increase MAXRING') + 30 format (/,' RINGS -- Too many 4-Membered Rings') call fatal end if iring4(1,nring4) = ia @@ -165,8 +162,7 @@ subroutine rings nring5 = nring5 + 1 if (nring5 .gt. maxring) then write (iout,50) - 50 format (/,' RINGS -- Too many 5-Membered Rings;', - & ' Increase MAXRING') + 50 format (/,' RINGS -- Too many 5-Membered Rings') call fatal end if iring5(1,nring5) = ia @@ -226,8 +222,8 @@ subroutine rings nring6 = nring6 + 1 if (nring6 .gt. maxring) then write (iout,70) - 70 format (/,' RINGS -- Too many 6-Membered', - & ' Rings; Increase MAXRING') + 70 format (/,' RINGS -- Too many', + & ' 6-Membered Rings') call fatal end if iring6(1,nring6) = ia diff --git a/source/scales.f b/source/scales.f index 7d4e16be6..1c8e0f2c4 100644 --- a/source/scales.f +++ b/source/scales.f @@ -19,7 +19,7 @@ module scales use sizes implicit none - real*8 scale(3*maxatm) + real*8 scale(maxvar) logical set_scale save end diff --git a/source/shakeup.f b/source/shakeup.f index bf813ddad..04b340136 100644 --- a/source/shakeup.f +++ b/source/shakeup.f @@ -385,6 +385,7 @@ subroutine shakeup c c subroutine chkangle (ia,ib,ic) + use sizes use couple use freeze use ring diff --git a/source/sizes.f b/source/sizes.f index 6a59ee35c..560638b1e 100644 --- a/source/sizes.f +++ b/source/sizes.f @@ -26,8 +26,12 @@ c maxref stored reference molecular systems c maxtyp force field atom type definitions c maxclass force field atom class definitions +c maxprm lines in the parameter file +c maxkey lines in the keyword file c maxrot bonds for torsional rotation +c maxvar optimization variables (vector storage) c maxopt optimization variables (matrix storage) +c maxlight sites for method of lights neighbors c maxvlst neighbors in van der Waals pair list c maxelst neighbors in electrostatics pair list c maxulst neighbors in dipole preconditioner list @@ -36,8 +40,16 @@ c maxvib vibrational frequencies c maxgeo distance geometry points c maxcell unit cells in replicated crystal +c maxring 3-, 4-, or 5-membered rings c maxbio biopolymer atom definitions c maxres residues in the macromolecule +c maxele elements in periodic table +c maxamino amino acid residue types +c maxnuc nucleic acid residue types +c maxbnd covalent bonds in molecular system +c maxang bond angles in molecular system +c maxtors torsional angles in molecular system +c maxbitor bitorsions in molecular system c c module sizes @@ -45,20 +57,30 @@ module sizes integer maxatm,maxval integer maxgrp,maxref integer maxtyp,maxclass - integer maxrot,maxopt + integer maxprm,maxkey + integer maxrot,maxvar + integer maxopt,maxlight integer maxvlst,maxelst integer maxulst,maxfft integer maxfix,maxvib integer maxgeo,maxcell - integer maxbio,maxres + integer maxring,maxbio + integer maxres,maxele + integer maxamino,maxnuc + integer maxbnd,maxang + integer maxtors,maxbitor parameter (maxatm=1000000) parameter (maxval=8) parameter (maxgrp=1000) parameter (maxref=10) parameter (maxtyp=5000) parameter (maxclass=1000) + parameter (maxprm=25000) + parameter (maxkey=25000) parameter (maxrot=1000) + parameter (maxvar=3*maxatm) parameter (maxopt=1000) + parameter (maxlight=8*maxatm) parameter (maxvlst=1800) parameter (maxelst=1200) parameter (maxulst=100) @@ -67,7 +89,15 @@ module sizes parameter (maxvib=1000) parameter (maxgeo=2500) parameter (maxcell=10000) + parameter (maxring=10000) parameter (maxbio=10000) parameter (maxres=10000) + parameter (maxele=112) + parameter (maxamino=38) + parameter (maxnuc=12) + parameter (maxbnd=2*maxatm) + parameter (maxang=4*maxatm) + parameter (maxtors=6*maxatm) + parameter (maxbitor=8*maxatm) save end diff --git a/source/solute.f b/source/solute.f index fa0570154..296661a48 100644 --- a/source/solute.f +++ b/source/solute.f @@ -13,17 +13,17 @@ c c c doffset dielectric offset to continuum solvation atomic radii -c p1 single-atom scale factor for analytical Still radii -c p2 1-2 interaction scale factor for analytical Still radii -c p3 1-3 interaction scale factor for analytical Still radii -c p4 nonbonded scale factor for analytical Still radii -c p5 soft cutoff parameter for analytical Still radii c rsolv atomic radius of each atom for continuum solvation c asolv atomic surface area solvation parameters c rborn Born radius of each atom for GB/SA solvation c drb solvation derivatives with respect to Born radii c drbp GK polarization derivatives with respect to Born radii c drobc chain rule term for Onufriev-Bashford-Case radii +c p1 single-atom scale factor for analytical Still radii +c p2 1-2 interaction scale factor for analytical Still radii +c p3 1-3 interaction scale factor for analytical Still radii +c p4 nonbonded scale factor for analytical Still radii +c p5 soft cutoff parameter for analytical Still radii c gpol polarization self-energy values for each atom c shct overlap scale factors for Hawkins-Cramer-Truhlar radii c aobc alpha values for Onufriev-Bashford-Case radii @@ -38,24 +38,25 @@ c c module solute + use sizes implicit none real*8 doffset + real*8 rsolv(maxatm) + real*8 asolv(maxatm) + real*8 rborn(maxatm) + real*8 drb(maxatm) + real*8 drbp(maxatm) + real*8 drobc(maxatm) real*8 p1,p2,p3,p4,p5 - real*8, allocatable :: rsolv(:) - real*8, allocatable :: asolv(:) - real*8, allocatable :: rborn(:) - real*8, allocatable :: drb(:) - real*8, allocatable :: drbp(:) - real*8, allocatable :: drobc(:) - real*8, allocatable :: gpol(:) - real*8, allocatable :: shct(:) - real*8, allocatable :: aobc(:) - real*8, allocatable :: bobc(:) - real*8, allocatable :: gobc(:) - real*8, allocatable :: vsolv(:) - real*8, allocatable :: wace(:,:) - real*8, allocatable :: s2ace(:,:) - real*8, allocatable :: uace(:,:) + real*8 gpol(maxatm) + real*8 shct(maxatm) + real*8 aobc(maxatm) + real*8 bobc(maxatm) + real*8 gobc(maxatm) + real*8 vsolv(maxatm) + real*8 wace(maxclass,maxclass) + real*8 s2ace(maxclass,maxclass) + real*8 uace(maxclass,maxclass) character*8 solvtyp character*8 borntyp save diff --git a/source/stochisok.f b/source/stochisok.f new file mode 100644 index 000000000..346523c56 --- /dev/null +++ b/source/stochisok.f @@ -0,0 +1,486 @@ +c +c +c ################################################### +c ## COPYRIGHT (C) 2011 by Jay William Ponder ## +c ## All Rights Reserved ## +c ################################################### +c +c ############################################################# +c ## ## +c ## subroutine stochisok ## +c ## ## +c ############################################################# +c +c +c + subroutine stochisok (istep,dt) + use sizes + use atomid + use freeze + use moldyn + use units + use virial + use atoms + use bath + use bound + use inform + use iounit + use keys + use mdstuf + use potent + use solute + use stodyn + use usage + implicit none + integer i,j,k + integer istep + real*8 dt,dt_2 + real*8 dta,dta_2 + real*8 epot,etot + real*8 eksum,eps + real*8 temp,pres + real*8 ealt,dalt + real*8 kBT,LkT,tempstor,temp2,len_fac + real*8 ekin(3,3) + real*8 stress(3,3) + real*8 viralt(3,3) + real*8, allocatable :: xold(:) + real*8, allocatable :: yold(:) + real*8, allocatable :: zold(:) + real*8, allocatable :: derivs(:,:) + + + +c +c +c set some values for the dynamics integration +c + dt_2 = 0.5d0 * dt + dta = dt / (dble(nrespa)) + dta_2 = 0.5d0 * dta + + kBT = boltzmann*kelvin + LkT = len_nhc*kBT + len_fac=dble(len_nhc)/dble(len_nhc+1) +c +c +c perform dynamic allocation of some local arrays +c + allocate (derivs(3,n)) + allocate (xold(n)) + allocate (yold(n)) + allocate (zold(n)) + +c +c +c initialize virial from fast-evolving potential energy terms +c + do i = 1, 3 + do j = 1, 3 + viralt(j,i) = 0.0d0 + end do + end do + +c +c +c begin inner RESPA loop +c + + do k = 1, nrespa + do i = 1, n + if (use(i)) then + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + end if + end do + +c +c +c nose-like isokinetic operator +c + call applyisok (dta) +c +c +c stochastic noise operator +c + call v2random (dta) +c +c +c force-like isokinetic operator +c + + call applyisok2 (dta) +c +c Evolve positions +c + if (use_rattle) then + do i = 1, n + xold(i) = x(i) + yold(i) = y(i) + zold(i) = z(i) + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end do + call rattle (dta,xold,yold,zold) + else + do i = 1, n + if (use(i)) then + x(i) = x(i) + v(1,i)*dta + y(i) = y(i) + v(2,i)*dta + z(i) = z(i) + v(3,i)*dta + end if + end do + end if +c +c get the potential energy and atomic forces +c + if (k.eq.nrespa) then + call gradfast (ealt,derivs) + + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call gradslow (epot,derivs) + do i = 1, n + do j = 1, 3 + a(j,i) = a(j,i) - ( nrespa * convert + & * derivs(j,i) / mass(i) ) + end do + end do + end if + + else + call gradfast (ealt,derivs) + if (use_rattle) then + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + call rattle2 (dta) + else + do i = 1, n + if (use(i)) then + do j = 1, 3 + a(j,i) = -convert * derivs(j,i) / mass(i) + end do + end if + end do + end if + end if + +c +c +c force-like isokinetic operator +c + call applyisok2 (dta) +c +c +c nose-like isokinetic operator +c + call applyisok (dta) + + end do + +c +c find the constraint-corrected full-step velocities +c + if (use_rattle) call rattle2 (dt) +c +c total potential and virial from sum of fast and slow parts +c + epot = epot + ealt + do i = 1, 3 + do j = 1, 3 + vir(j,i) = vir(j,i) + viralt(j,i) + end do + end do + + +c +c perform deallocation of some local arrays +c + + deallocate (derivs) + deallocate (xold) + deallocate (yold) + deallocate (zold) +c +c total energy is sum of kinetic and potential energies +c + etot = eksum + epot + + call kinetic (eksum,ekin) + call pressure (dt,epot,ekin,temp,pres,stress) + temp = 2.0d0 * eksum / (dble(nfree) * gasconst) + + +c tempstor=0.0d0 +c do i = 1, n +c do j = 1, 3 +c tempstor=tempstor+(mass(i)*v(j,i)*v(j,i)) +c do k = 1, len_nhc +c tempstor=tempstor+len_fac*q_iso1(k,j,i) +c & *v_iso1(k,j,i)*v_iso1(k,j,i) +c end do +c end do +c end do +c +c print (actual/target) isokinetic constraint +c +c write (iout,10) tempstor/(3.0d0*n*LkT) +c 10 format (3f16.14) + + + +c +c compute statistics and save trajectory for this step +c + call mdstat (istep,dt,etot,epot,eksum,temp,pres) + call mdsave (istep,dt,epot) + call mdrest (istep) + return + end + + +c +c ################################################################## +c ## ## +c ## subroutine applyisok ## +c ## ## +c ################################################################## +c +c +c "applyisok" represents the Liouville operator iL_N in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine applyisok (dt) + use sizes + use atomid + use atoms + use bath + use freeze + use moldyn + use units + use usage + use virial + implicit none + integer nc,ns,i,j,k,iatom,inhc,ixyz + real*8 dt,dtc,dts + real*8 dt2,dt4,dt8 + real*8 tempstor,kBT,LkT + real*8 w(3) + real*8 len_fac + + len_fac = dble(len_nhc)/dble(len_nhc+1) + + kBT = boltzmann*kelvin + LkT = (dble(len_nhc))*kBT + + nc = 5 + ns = 3 + dtc = dt / dble(nc) + w(1) = 1.0d0 / (2.0d0-2.0d0**(1.0d0/3.0d0)) + w(2) = 1.0d0 - 2.0d0*w(1) + w(3) = w(1) + +c +c use multiple time steps and Suzuki-Yoshida decomposition +c + do k = 1, nc + do j = 1, ns + + dts = w(j) * dtc + dt2 = 0.5d0 * dts + dt4 = 0.25d0 * dts + dt8 = 0.125d0 * dts + + do iatom = 1, n + do inhc = 1, len_nhc + do ixyz = 1, 3 + tempstor = q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc,ixyz,iatom)*v_iso1(inhc,ixyz,iatom) + & - kBT + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) + & + dt4*tempstor/q_iso2(inhc,ixyz,iatom) + end do + end do + end do + + do iatom = 1, n + do ixyz = 1, 3 + tempstor = mass(iatom)*v(ixyz,iatom)*v(ixyz,iatom) + do inhc = 1, len_nhc + tempstor = tempstor + (len_fac * + & q_iso1(inhc,ixyz,iatom) * + & v_iso1(inhc,ixyz,iatom) * + & v_iso1(inhc,ixyz,iatom) * + & exp(-2.0 * dt2 * v_iso2(inhc,ixyz,iatom))) + end do + tempstor = sqrt(LkT/tempstor) + v(ixyz,iatom) = v(ixyz,iatom) * tempstor + do inhc = 1, len_nhc + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom) * + & tempstor*exp(-dt2*v_iso2(inhc,ixyz,iatom)) + end do + end do + end do + + do iatom = 1, n + do inhc = 1, len_nhc + do ixyz = 1, 3 + tempstor = q_iso1(inhc,ixyz,iatom) + & *v_iso1(inhc,ixyz,iatom)*v_iso1(inhc,ixyz,iatom) + & - kBT + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) + & + dt4*tempstor/q_iso2(inhc,ixyz,iatom) + end do + end do + end do + + end do + end do +c do i = 1, 1 +c do j = 1, 1 +c write (*,10) v(j,i) +c do k = 1, len_nhc +c write (*,10) v_iso1(k,j,i) +c end do +c end do +c end do +10 format (3f16.6) + return + end +c +c ################################################################## +c ## ## +c ## subroutine applyisok2 ## +c ## ## +c ################################################################## +c +c +c "applyisok2" represents the Liuville operator iL_v in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine applyisok2 (dt) + use sizes + use atomid + use atoms + use bath + use freeze + use moldyn + use units + use usage + use virial + implicit none + real*8 kBT,LkT,isoA,isoB + real*8 dt,dt_2,rootB,arg + real*8 iso_s,iso_sdot + real*8 epot,etot + integer iatom,inhc,ixyz,i,j,k + + dt_2 = 0.5d0*dt + kBT = boltzmann*kelvin + LkT = (dble(len_nhc))*kBT + + + do iatom = 1, n + do ixyz = 1, 3 + isoA = mass(iatom)*a(ixyz,iatom)*v(ixyz,iatom)/LkT + isoB = mass(iatom)*a(ixyz,iatom)*a(ixyz,iatom)/LkT + rootB = sqrt(isoB) + arg = dt_2*rootB +c write (*,10) isoA +c write (*,10) isoB +c write (*,10) LkT + if (arg > 0.00001) then + iso_s = (1.0/rootB)*sinh(arg) + + & (isoA/isoB)*(cosh(arg)-1.0) + iso_sdot = cosh(arg) + (isoA/rootB)*sinh(arg) + else + iso_s = ((((isoB*isoA/24.0)*dt_2 +isoB/6.0)*dt_2 + + & 0.5*isoA)*dt_2 + 1.0)*dt_2 + iso_sdot = (((isoB*isoA/6.0)*dt_2 +0.5*isoB)*dt_2 + + & isoA)*dt_2 + 1.0 + end if + v(ixyz,iatom) = (v(ixyz,iatom) + (a(ixyz,iatom))*iso_s) + & /iso_sdot + do inhc = 1, len_nhc + v_iso1(inhc,ixyz,iatom)=v_iso1(inhc,ixyz,iatom)/iso_sdot + end do + end do + end do + + + return + end +c +c ################################################################## +c ## ## +c ## subroutine v2random ## +c ## ## +c ################################################################## +c +c +c "v2random" represents the Liuville operator iL_OU in Leimkuhler, Margul, and Tuckerman (2013) +c +c + subroutine v2random (dt) + use sizes + use atomid + use atoms + use bath + use freeze + use moldyn + use units + use usage + use virial + implicit none + real*8 dt,dt_2,kBT + real*8 gamma,stoch_e,sigma + real*8 normal + integer iatom,inhc,ixyz + + kBT = boltzmann*kelvin + gamma = 1.0d0 / (200.0d0*dt) + sigma = sqrt(kBT* (1.0d0-exp(-2.0d0*gamma*dt) )/q_iso2(1,1,1) ) + stoch_e = exp(-gamma*dt) + + do iatom = 1, n + do ixyz = 1, 3 + do inhc = 1, len_nhc + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) * + & stoch_e + v_iso2(inhc,ixyz,iatom) = v_iso2(inhc,ixyz,iatom) + + & (sigma * normal () ) + end do + end do + end do + + return + end \ No newline at end of file diff --git a/source/testgrad.f b/source/testgrad.f index 15c16fc1e..8539e7372 100644 --- a/source/testgrad.f +++ b/source/testgrad.f @@ -127,7 +127,7 @@ program testgrad 50 continue if (query) then write (iout,60) eps0 - 60 format (/,' Enter Finite Difference Stepsize [',d8.1, + 60 format (/,' Enter a Numerical Stepsize [',d7.1, & ' Ang] : ',$) read (input,70,err=50) eps 70 format (f20.0) diff --git a/source/testhess.f b/source/testhess.f index f950073ba..f8c24f82b 100644 --- a/source/testhess.f +++ b/source/testhess.f @@ -133,7 +133,7 @@ program testhess 70 continue if (query) then write (iout,80) eps0 - 80 format (/,' Enter Finite Difference Stepsize [',d8.1, + 80 format (/,' Enter a Numerical Stepsize [',d7.1, & ' Ang] : ',$) read (input,90,err=70) eps 90 format (f20.0) diff --git a/source/testpair.f b/source/testpair.f index 01e790b6a..414cafb66 100644 --- a/source/testpair.f +++ b/source/testpair.f @@ -216,12 +216,12 @@ program testpair c c perform dynamic allocation of some global arrays c - allocate (dev(3,n)) - allocate (dec(3,n)) - allocate (decd(3,n)) - allocate (ded(3,n)) - allocate (dem(3,n)) - allocate (dep(3,n)) + allocate (dev(3,maxatm)) + allocate (dec(3,maxatm)) + allocate (decd(3,maxatm)) + allocate (ded(3,maxatm)) + allocate (dem(3,maxatm)) + allocate (dep(3,maxatm)) c c get the timing for energy terms via double nested loop c diff --git a/source/testrot.f b/source/testrot.f index 76bc9d409..37d6a3308 100644 --- a/source/testrot.f +++ b/source/testrot.f @@ -63,7 +63,7 @@ program testrot 10 continue if (query) then write (iout,20) delta0 - 20 format (/,' Enter Finite Difference Stepsize [',d8.1, + 20 format (/,' Enter Numerical Derivative Stepsize [',d7.1, & ' Deg] : ',$) read (input,30,err=10) delta 30 format (f20.0) diff --git a/source/tncg.f b/source/tncg.f index f63c59130..6dfb31394 100644 --- a/source/tncg.f +++ b/source/tncg.f @@ -144,6 +144,12 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, c c check number of variables and get type of optimization c + if (nvar .gt. maxvar) then + write (iout,10) + 10 format (/,' TNCG -- Too many Parameters,', + & ' Increase the Value of MAXVAR') + return + end if rms = sqrt(dble(nvar)) if (coordtype .eq. 'CARTESIAN') then rms = rms / sqrt(3.0d0) @@ -186,29 +192,29 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, call upcase (keyword) string = record(next:120) if (keyword(1:7) .eq. 'FCTMIN ') then - read (string,*,err=10,end=10) fctmin + read (string,*,err=20,end=20) fctmin else if (keyword(1:8) .eq. 'MAXITER ') then - read (string,*,err=10,end=10) maxiter + read (string,*,err=20,end=20) maxiter else if (keyword(1:9) .eq. 'NEXTITER ') then - read (string,*,err=10,end=10) nextiter + read (string,*,err=20,end=20) nextiter else if (keyword(1:8) .eq. 'NEWHESS ') then - read (string,*,err=10,end=10) newhess + read (string,*,err=20,end=20) newhess else if (keyword(1:12) .eq. 'SADDLEPOINT ') then negtest = .false. else if (keyword(1:8) .eq. 'STEPMIN ') then - read (string,*,err=10,end=10) stpmin + read (string,*,err=20,end=20) stpmin else if (keyword(1:8) .eq. 'STEPMAX ') then - read (string,*,err=10,end=10) stpmax + read (string,*,err=20,end=20) stpmax else if (keyword(1:6) .eq. 'CAPPA ') then - read (string,*,err=10,end=10) cappa + read (string,*,err=20,end=20) cappa else if (keyword(1:9) .eq. 'SLOPEMAX ') then - read (string,*,err=10,end=10) slpmax + read (string,*,err=20,end=20) slpmax else if (keyword(1:7) .eq. 'ANGMAX ') then - read (string,*,err=10,end=10) angmax + read (string,*,err=20,end=20) angmax else if (keyword(1:7) .eq. 'INTMAX ') then - read (string,*,err=10,end=10) intmax + read (string,*,err=20,end=20) intmax end if - 10 continue + 20 continue end do c c initialize the function call and iteration counters @@ -222,27 +228,27 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, c if (iprint .gt. 0) then if (mode .eq. 'NEWTON') then - write (iout,20) - 20 format (/,' Full-Newton Conjugate-Gradient', + write (iout,30) + 30 format (/,' Full-Newton Conjugate-Gradient', & ' Optimization :') else if (mode .eq. 'TNCG') then - write (iout,30) - 30 format (/,' Truncated-Newton Conjugate-Gradient', + write (iout,40) + 40 format (/,' Truncated-Newton Conjugate-Gradient', & ' Optimization :') else if (mode .eq. 'DTNCG') then - write (iout,40) - 40 format (/,' Finite-Difference Truncated-Newton', + write (iout,50) + 50 format (/,' Finite-Difference Truncated-Newton', & ' Conjugate-Gradient Optimization :') else if (mode .eq. 'AUTO') then - write (iout,50) - 50 format (/,' Variable-Mode Truncated-Newton', + write (iout,60) + 60 format (/,' Variable-Mode Truncated-Newton', & ' Conjugate-Gradient Optimization :') end if - write (iout,60) mode,method,grdmin - 60 format (/,' Algorithm : ',a6,5x,'Preconditioning : ',a6,5x, + write (iout,70) mode,method,grdmin + 70 format (/,' Algorithm : ',a6,5x,'Preconditioning : ',a6,5x, & ' RMS Grad :',d9.2) - write (iout,70) - 70 format (/,' TN Iter F Value G RMS F Move ', + write (iout,80) + 80 format (/,' TN Iter F Value G RMS F Move ', & ' X Move CG Iter Solve FG Call',/) end if c @@ -277,11 +283,11 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, c if (iprint .gt. 0) then if (f.lt.1.0d7 .and. f.gt.-1.0d6 .and. g_rms.lt.1.0d6) then - write (iout,80) iter_tn,f,g_rms,fg_call - 80 format (i6,f13.4,f12.4,41x,i7) - else write (iout,90) iter_tn,f,g_rms,fg_call - 90 format (i6,d13.4,d12.4,41x,i7) + 90 format (i6,f13.4,f12.4,41x,i7) + else + write (iout,100) iter_tn,f,g_rms,fg_call + 100 format (i6,d13.4,d12.4,41x,i7) end if end if c @@ -295,22 +301,22 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, done = .true. minimum = f if (iprint .gt. 0) then - write (iout,100) - 100 format (/,' TNCG -- Normal Termination due to SmallGrad') + write (iout,110) + 110 format (/,' TNCG -- Normal Termination due to SmallGrad') end if else if (f .le. fctmin) then done = .true. minimum = f if (iprint .gt. 0) then - write (iout,110) - 110 format (/,' TNCG -- Normal Termination due to SmallFct') + write (iout,120) + 120 format (/,' TNCG -- Normal Termination due to SmallFct') end if else if (iter_tn .ge. maxiter) then done = .true. minimum = f if (iprint .gt. 0) then - write (iout,120) - 120 format (/,' TNCG -- Incomplete Convergence', + write (iout,130) + 130 format (/,' TNCG -- Incomplete Convergence', & ' due to IterLimit') end if end if @@ -445,13 +451,13 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, if (done .or. mod(iter_tn,iprint).eq.0) then if (f.lt.1.0d7 .and. f.gt.-1.0d6 .and. & g_rms.lt.1.0d6 .and. f_move.lt.1.0d5) then - write (iout,130) iter_tn,f,g_rms,f_move,x_move, + write (iout,140) iter_tn,f,g_rms,f_move,x_move, & iter_cg,info_solve,fg_call - 130 format (i6,f13.4,f12.4,f11.4,f10.4,i8,3x,a9,i7) + 140 format (i6,f13.4,f12.4,f11.4,f10.4,i8,3x,a9,i7) else - write (iout,140) iter_tn,f,g_rms,f_move,x_move, + write (iout,150) iter_tn,f,g_rms,f_move,x_move, & iter_cg,info_solve,fg_call - 140 format (i6,d13.4,d12.4,d11.4,f10.4,i8,3x,a9,i7) + 150 format (i6,d13.4,d12.4,d11.4,f10.4,i8,3x,a9,i7) end if end if end if @@ -470,11 +476,11 @@ subroutine tncg (mode,method,nvar,x0,minimum,grdmin, minimum = f if (iprint .gt. 0) then if (g_rms.le.grdmin .or. f.le.fctmin) then - write (iout,150) status - 150 format (/,' TNCG -- Normal Termination due to ',a9) - else write (iout,160) status - 160 format (/,' TNCG -- Incomplete Convergence', + 160 format (/,' TNCG -- Normal Termination due to ',a9) + else + write (iout,170) status + 170 format (/,' TNCG -- Incomplete Convergence', & ' due to ',a9) end if end if diff --git a/source/torsions.f b/source/torsions.f index 923ce5cd9..550b62832 100644 --- a/source/torsions.f +++ b/source/torsions.f @@ -18,7 +18,6 @@ c subroutine torsions use sizes - use atoms use bndstr use couple use iounit @@ -26,12 +25,10 @@ subroutine torsions implicit none integer i,j,k integer ia,ib,ic,id - integer maxtors c c c perform dynamic allocation of some global arrays c - maxtors = 6 * n if (.not. allocated(itors)) allocate (itors(4,maxtors)) c c loop over all bonds, storing the atoms in each torsion diff --git a/source/vibbig.f b/source/vibbig.f index eb5b51f3c..205c6f186 100644 --- a/source/vibbig.f +++ b/source/vibbig.f @@ -2071,9 +2071,9 @@ subroutine hessblk (amass,k0,i1,i2,vector) c if (first) then first = .false. - if (.not. allocated(hessx)) allocate (hessx(3,n)) - if (.not. allocated(hessy)) allocate (hessy(3,n)) - if (.not. allocated(hessz)) allocate (hessz(3,n)) + if (.not. allocated(hessx)) allocate (hessx(3,maxatm)) + if (.not. allocated(hessy)) allocate (hessy(3,maxatm)) + if (.not. allocated(hessz)) allocate (hessz(3,maxatm)) end if c c zero out the Hessian elements for the current atom diff --git a/source/warp.f b/source/warp.f index 5607f5ad2..9edd7263f 100644 --- a/source/warp.f +++ b/source/warp.f @@ -13,10 +13,10 @@ c c c deform value of the smoothing deformation parameter +c m2 second moment of the GDA gaussian for each atom c difft diffusion coefficient for torsional potential c diffv diffusion coefficient for van der Waals potential c diffc diffusion coefficient for charge-charge potential -c m2 second moment of the GDA gaussian for each atom c use_smooth flag to use a potential energy smoothing method c use_dem flag to use diffusion equation method potential c use_gda flag to use gaussian density annealing potential @@ -25,12 +25,13 @@ c c module warp + use sizes implicit none real*8 deform + real*8 m2(maxatm) real*8 difft real*8 diffv real*8 diffc - real*8, allocatable :: m2(:) logical use_smooth logical use_dem logical use_gda