source/unres/src_MD/md-diff/np/nhcint.f

   1 c-----------------------------------------------------------------
   2       SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
   3       implicit real*8 (a-h,o-z)
   4       include 'DIMENSIONS'
   5       include 'COMMON.MD'
   6       double precision akin,gnkt,dt,aa,gkt,scale
   7       double precision wdti(maxyosh),wdti2(maxyosh),
   8      &                 wdti4(maxyosh),wdti8(maxyosh)
   9       integer i,iresn,iyosh,inos,nnos1
  10
  11       dt=d_time
  12       nnos1=nnos+1
  13       GKT = Rb*t_bath
  14       GNKT = dimen*GKT
  15       akin=akin*2
  16
  17
  18 C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
  19 C INTEGRATION FROM t=0 TO t=DT/2
  20 C GET THE TOTAL KINETIC ENERGY
  21       SCALE = 1.D0
  22 c      CALL GETKINP(MASS,VX,VY,VZ,AKIN)
  23 C UPDATE THE FORCES
  24       GLOGS(1) = (AKIN - GNKT)/QMASS(1)
  25 C START THE MULTIPLE TIME STEP PROCEDURE
  26       DO IRESN = 1,NRESN
  27        DO IYOSH = 1,NYOSH
  28 C UPDATE THE THERMOSTAT VELOCITIES
  29         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
  30         DO INOS = 1,NNOS-1
  31          AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
  32          VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
  33      &          + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
  34         ENDDO
  35 C UPDATE THE PARTICLE VELOCITIES
  36         AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
  37         SCALE = SCALE*AA
  38 C UPDATE THE FORCES
  39         GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
  40 C UPDATE THE THERMOSTAT POSITIONS
  41         DO INOS = 1,NNOS
  42          XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
  43         ENDDO
  44 C UPDATE THE THERMOSTAT VELOCITIES
  45         DO INOS = 1,NNOS-1
  46          AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
  47          VLOGS(INOS) = VLOGS(INOS)*AA*AA
  48      &      + WDTI4(IYOSH)*GLOGS(INOS)*AA
  49          GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
  50      &      -GKT)/QMASS(INOS+1)
  51         ENDDO
  52         VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
  53        ENDDO
  54       ENDDO
  55 C UPDATE THE PARTICLE VELOCITIES
  56 c outside of this subroutine
  57 c      DO I = 1,N
  58 c       VX(I) = VX(I)*SCALE
  59 c       VY(I) = VY(I)*SCALE
  60 c       VZ(I) = VZ(I)*SCALE
  61 c      ENDDO
  62       RETURN
  63       END