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