+++ /dev/null
-c-----------------------------------------------------------------
- SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.MD'
- double precision akin,gnkt,dt,aa,gkt,scale
- double precision wdti(maxyosh),wdti2(maxyosh),
- & wdti4(maxyosh),wdti8(maxyosh)
- integer i,iresn,iyosh,inos,nnos1
-
- dt=d_time
- nnos1=nnos+1
- GKT = Rb*t_bath
- GNKT = dimen*GKT
- akin=akin*2
-
-
-C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
-C INTEGRATION FROM t=0 TO t=DT/2
-C GET THE TOTAL KINETIC ENERGY
- SCALE = 1.D0
-c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
-C UPDATE THE FORCES
- GLOGS(1) = (AKIN - GNKT)/QMASS(1)
-C START THE MULTIPLE TIME STEP PROCEDURE
- DO IRESN = 1,NRESN
- DO IYOSH = 1,NYOSH
-C UPDATE THE THERMOSTAT VELOCITIES
- VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
- DO INOS = 1,NNOS-1
- AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
- VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
- & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
- ENDDO
-C UPDATE THE PARTICLE VELOCITIES
- AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
- SCALE = SCALE*AA
-C UPDATE THE FORCES
- GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
-C UPDATE THE THERMOSTAT POSITIONS
- DO INOS = 1,NNOS
- XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
- ENDDO
-C UPDATE THE THERMOSTAT VELOCITIES
- DO INOS = 1,NNOS-1
- AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
- VLOGS(INOS) = VLOGS(INOS)*AA*AA
- & + WDTI4(IYOSH)*GLOGS(INOS)*AA
- GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
- & -GKT)/QMASS(INOS+1)
- ENDDO
- VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
- ENDDO
- ENDDO
-C UPDATE THE PARTICLE VELOCITIES
-c outside of this subroutine
-c DO I = 1,N
-c VX(I) = VX(I)*SCALE
-c VY(I) = VY(I)*SCALE
-c VZ(I) = VZ(I)*SCALE
-c ENDDO
- RETURN
- END