1 c-----------------------------------------------------------------
2 SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
3 implicit real*8 (a-h,o-z)
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
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
22 c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
24 GLOGS(1) = (AKIN - GNKT)/QMASS(1)
25 C START THE MULTIPLE TIME STEP PROCEDURE
28 C UPDATE THE THERMOSTAT VELOCITIES
29 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
31 AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
32 VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
33 & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
35 C UPDATE THE PARTICLE VELOCITIES
36 AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
39 GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
40 C UPDATE THE THERMOSTAT POSITIONS
42 XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
44 C UPDATE THE THERMOSTAT VELOCITIES
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)
52 VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
55 C UPDATE THE PARTICLE VELOCITIES
56 c outside of this subroutine