C 22 Feb 06 - HL  - fix parallel print out
C 14 Nov 05 - HMN - NVT MD run with temperature reassignment
C 19 Sep 05 - HMN,IA - implement EFP minimum image convention MD
C
C*MODULE MDDRVR   *DECK MDX
      SUBROUTINE MDX
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00)
      PARAMETER (MXATM=2000)
      PARAMETER (MXFRG=1050)
      PARAMETER (MXFGPT=12000)
      PARAMETER (MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (NBIN=1000)
C
      LOGICAL READ,NVT,MBT,MBR,QRAND,EQUIL,PBCOUT,PROD
      LOGICAL PBCWRI
C
      DOUBLE PRECISION MD,MDTYP,MDINT
      INTEGER TYPMD, INTMD, RAT,RETEMP
C---      integer pevery, pt
C
      DIMENSION TVEL( 3, MXFRG )
      DIMENSION RMOM( 3, MXFRG ), RVEL( 3, MXFRG )
      DIMENSION QUAT(4, MXFRG)
      DIMENSION QUAT1D( 4, MXFRG ), QUAT2D( 4, MXFRG )
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /FUNCT / E, EG(3*MXATM)
      COMMON /GOFR  / GOFR(NBIN,3), DELR, NUMG, NUMSAM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /MDSIM2/ READ
      COMMON /MDSIM3/ NVT
      COMMON /MDSIM4/ NSTEPS
      COMMON /MDSIM5/ MBT,MBR,QRAND
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PBC/    PBCOUT,PBCWRI
C
C         set up namelist simulation
      PARAMETER (NNAM=24)
      DIMENSION QNAM( NNAM )
      DIMENSION KQNAM( NNAM )
      DATA MD/8HMD      /
      DATA QNAM /8HMDTYP   ,8HMDINT   ,8HBATHT   ,8HBOXL    ,
     *           8HDELR    ,8HDT      ,8HJEVERY  ,8HKEVERY  ,
     *           8HRETEMP  ,8HNSTEPS  ,8HNVT     ,
     *           8HREAD    ,8HMBT     ,8HMBR     ,8HQRAND   ,
     *           8HPROD    ,
     *           8HEQUIL   ,8HPBCOUT  ,8HTVEL    ,8HRMOM    ,
     *           8HRVEL    ,8HQUAT    ,8HQUAT1D  ,8HQUAT2D  /
      DATA KQNAM/5, 5, 3, 3, 3, 3, 1, 1, 1, 1, 0, 0, 0, 0, 0,
     *                   0, 0, 0, -3, -3, -3, -3, -3, -3/
C
      DATA BLANK      /8H        /
      DATA EFP /8HEFP     /
      DATA FROG/8HFROG    /
C
C*********************************************************
C** DRIVER FOR MD SIMULATIONS USING EFFECTIVE FRAGMENTS **
C**           JON SORENSON, 12/7/99        **
C**  (DRIVER CODE BORROWED HEAVILY FROM GRANT MERRILL)  **
C*********************************************************
C
C  Insert another parameter into input: septol = separation
C  tolerance between atoms of different molecules--used in
C  polarization and exchange repulsion calculations.
C
C   *** set default variable values for $md group ***
C
      MDTYP = BLANK
      MDINT = BLANK
      BATHT = 300.00D+00
      BOXL = 12.43D+00
      BOXLI = ONE / BOXL
      DELR = 0.02D+00
      DT = 1.0D-15
C
      JEVERY = 10
      KEVERY = 100
      RETEMP = 0
      NSTEPS = 10000
C       nmax = 5
C       nsqmax = nmax * nmax + 1
C       septol = 1.30d+00
      NVT    = .FALSE.
      READ = .FALSE.
      MBT = .FALSE.
      MBR = .FALSE.
      QRAND = .FALSE.
      EQUIL = .FALSE.
      PBCOUT = .FALSE.
      PBCWRI = .FALSE.
      PROD = .FALSE.
C
      CALL VCLR(TVEL,1,3*MXFRG)
      CALL VCLR(QUAT,1,4*MXFRG)
      CALL VCLR(RVEL,1,3*MXFRG)
      CALL VCLR(RMOM,1,3*MXFRG)
      CALL VCLR(QUAT1D,1,4*MXFRG)
      CALL VCLR(QUAT2D,1,4*MXFRG)
C
      DO I=1,NFRG
         QUAT(1,I) = 1.0D+00
      ENDDO
C
C       *** open external .irc file ***
C
      IPIRC = 4
      CALL SEQOPN( IPIRC, 'IRCDATA', 'NEW', .FALSE., 'FORMATTED')
C
C       *** read user-defined variables from $md group
C
C   Vectors must be treated differently
      KQNAM(19)= 100*MXFRG + 3
      KQNAM(20)= 100*MXFRG + 3
      KQNAM(21)= 100*MXFRG + 3
      KQNAM(22)= 100*MXFRG + 3
      KQNAM(23)= 100*MXFRG + 3
      KQNAM(24)= 100*MXFRG + 3
      IRET = 0
      CALL NAMEIO(IR,IRET,MD,NNAM,QNAM,KQNAM,
     *            MDTYP,MDINT,BATHT,BOXL,DELR,DT,
     *            JEVERY,KEVERY,RETEMP,NSTEPS,
     *            NVT,READ,MBT,MBR,QRAND,PROD,EQUIL,PBCOUT,
     *            TVEL,RMOM,RVEL,QUAT,QUAT1D,QUAT2D,
     *            0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *            0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      IF ( IRET .EQ. 1 .OR. IRET .EQ. 2 ) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR READING $MD GROUP'
         CALL ABRT
      ENDIF
C
      NERR = 0
      IF (MDTYP .EQ. BLANK) THEN
         WRITE(IW,4999)
         NERR = NERR + 1
      ENDIF
C
      CALL MDCONV
C
      BOXL = BOXL / AU2ANG
      BOXLI = ONE / BOXL
C
C        *** check some variable values ***
C
        IF ( DT .GT. 5.0D-15 ) THEN
          WRITE(IW,5000)
          NERR = NERR + 1
        ENDIF
C
        IF ( NERR .GT. 0 ) THEN
          WRITE(IW,5050) NERR
          WRITE(IW,5060)
          CALL ABRT
        ENDIF
C
C
C        *** output $md variable values ***
C
        IF (MASWRK) THEN
          WRITE(IW,6000) MDTYP, MDINT, BATHT, BOXL*AU2ANG,
     *                   DELR, DT, JEVERY, KEVERY,
     *                   RETEMP, NSTEPS, NVT, READ,
     *                   MBT, MBR, QRAND, PROD, EQUIL, PBCOUT
          IF(MDTYP.EQ.EFP) THEN
            WRITE(IW,6010) NFRG
          ELSE
            WRITE(IW,6011) NAT
          ENDIF
        ENDIF
C
       IF(MDTYP.EQ.EFP) THEN
        TYPMD = 1
       ELSE
        TYPMD = 2
       ENDIF
C
       IF(MDINT.EQ.FROG) THEN
        INTMD = 1
       ELSE
        INTMD = 2
       ENDIF
C
C        *** calculate initial energy and forces ***
C
        IF (TYPMD.EQ.1) THEN
             CALL EFCM
             CALL TIMIT(1)
             CALL ENEDEF
             CALL TIMIT(1)
        ELSE
           CALL GRADX
        ENDIF
C
C        *** initialize various quantities, convert to quaternions ***
C
C--For EFPs:
        IF (TYPMD.EQ.1) THEN
               CALL MODSET
C
C--If reading:
           IF(READ) THEN
              FACTOR = 1.0D+00 / 1.0D-12
              DO I = 1,NFRG
                VX(I) = TVEL(1,I)*FACTOR
                VY(I) = TVEL(2,I)*FACTOR
                VZ(I) = TVEL(3,I)*FACTOR
                QW(I) = QUAT(1,I)
                QX(I) = QUAT(2,I)
                QY(I) = QUAT(3,I)
                QZ(I) = QUAT(4,I)
              ENDDO
C
C--For frog, setup ang. mom., for verl, set w and qdots..
              IF(INTMD.EQ.1) THEN
                DO I=1,NFRG
                 OX(I) = RMOM(1,I)*FACTOR
                 OY(I) = RMOM(2,I)*FACTOR
                 OZ(I) = RMOM(3,I)*FACTOR
                ENDDO
              ELSE
                DO I=1,NFRG
                 WX(I) = RVEL(1,I)*FACTOR
                 WY(I) = RVEL(1,I)*FACTOR
                 WZ(I) = RVEL(1,I)*FACTOR
                 QW1(I) = QUAT1D(1,I)
                 QX1(I) = QUAT1D(2,I)
                 QY1(I) = QUAT1D(3,I)
                 QZ1(I) = QUAT1D(4,I)
                 QW2(I) = QUAT2D(1,I)
                 QX2(I) = QUAT2D(2,I)
                 QY2(I) = QUAT2D(3,I)
                 QZ2(I) = QUAT2D(4,I)
                ENDDO
              ENDIF
           ELSE
C
C--No reading:
Cc
C--MB trans. vel...
              IF(MBT) THEN
               CALL NORMT(BATHT)
C               call gaussv(batht)
C---            call pchk(ptol)
              ELSE
               CALL CLRTVE
              ENDIF
C
C--MB rot. vel...
              IF(MBR) THEN
               CALL NORMW(BATHT)
C               call gaussw(batht)
              ELSE
               CALL CLRVEL
              ENDIF
C
C--Random quaternions...
              IF(QRAND) THEN
                CALL RANDQ
              ELSE
                CALL EU2QUT
              ENDIF
C
              IF(INTMD.EQ.1) THEN
               CALL CLRMOM
              ELSE
C               call clrqtd
               CALL INQDOT
              ENDIF
C
           ENDIF
        ELSE
           IF(READ) THEN
              FACTOR = 1.0D+00 / 1.0D-12
              DO I = 1,NAT
                VX(I) = TVEL(1,I)*FACTOR
                VY(I) = TVEL(2,I)*FACTOR
                VZ(I) = TVEL(3,I)*FACTOR
              ENDDO
            ELSE
              CALL CLRMOM
              CALL CLRTVE
            ENDIF
        ENDIF
        IF(PROD) CALL INITGR
        IF (MASWRK) THEN
          WRITE(IW,6020) E
          WRITE(IW,6030) E*AU2KCAL
C
C         write initial coords and velocities
          IF (TYPMD.EQ.1) THEN
             CALL MDWREF(0,INTMD)
             CALL MDRPREF(0)
          ELSE
             CALL MDWRAB(0)
             CALL MDRPRAB(0)
          ENDIF
          CALL TIMIT(1)
        ENDIF
C
C        *** start main simulation loop ***
        JT = 1
        KT = 1
        RAT = 1
        EKINT=ZERO
        EKINR=ZERO
C
C  Translational and Rotational KE at t
      DO I=1,NFRG
      EKINT = EKINT + TOTMAS(ISET(I))*(VX(I)**2 + VY(I)**2 + VZ(I)**2)
      EKINR = EKINR + AMOM(1,ISET(I)) * WX(I)**2 +
     *       AMOM(2,ISET(I)) * WY(I)**2 + AMOM(3,ISET(I)) * WZ(I)**2
        ENDDO
C
C          ekint = half * ek2kcal * totmas * ekint
           EKINT = HALF * EK2KCAL * EKINT
          EKINR = HALF * EK2KCAL * EKINR
C
C          write(6,*) 'initial KE: ekint,ekinr,ekin = ',
C     *        ekint,ekinr,ekint+ekinr
C          write(6,*) 'initial temp: ',tempmd(ekint+ekinr)
C
        DO I = 1, NSTEPS
          IF(PBCOUT.AND.I.EQ.NSTEPS) PBCWRI=.TRUE.
C
C        *** Seperate treatment for leapfrog and vverlet ***
C   Only calculate energy/grad outside if doing leapfrog;
C   energy/grad calculation done within integration for
C   vverlet...
C   Integrate equations of motion.
C
          IF(INTMD.EQ.1) THEN
           IF (TYPMD.EQ.1) THEN
                CALL ENEDEF
                CALL MDMOVE(INTMD)
           ELSE
            CALL ENEDAB
            CALL MDMOVA(EG,NAT,C)
           ENDIF
          ENDIF
C   Velocity verlet--need to change torque to torque/inertia for
C   first time through...
C
          IF(INTMD.EQ.2) THEN
           IF(I.EQ.1) THEN
            CALL TORQMO
           ENDIF
           CALL MOVEVV(INTMD)
          ENDIF
C
C       *** every retemp timesteps:                 ***
C       ***    if a nvt run, reassign temps from MB ***
C
            IF (NVT) THEN
             IF (RAT.EQ.RETEMP) THEN
              CALL NORMT(BATHT)
              CALL NORMW(BATHT)
              CALL EU2QUT
              CALL INQDOT
              RAT=0
              IF(MASWRK) WRITE(6,*) 'RESETTING VELOCITIES!!!!!!!!!'
             ENDIF
             RAT=RAT+1
            ENDIF
C
            IF (EQUIL) THEN
             IF (RAT.EQ.RETEMP) THEN
              CALL NORMT(BATHT)
              CALL NORMW(BATHT)
C              call eu2qut
C              call inqdot
              RAT=0
             ENDIF
             RAT=RAT+1
            ENDIF
C
C
C        *** every jevery timesteps:            ***
C        ***     report simulation quantities   ***
C        ***     collect g(r)s                  ***
C
C      scalt = sqrt(batht/(2.0d+00*ekint/(3.0d+00*nfrg)/1.987216d-03))
C      scalr = sqrt(batht/(2.0d+00*ekinr/(3.0d+00*nfrg)/1.987216d-03))
C
            IF ( JT .EQ. JEVERY ) THEN
             IF(PROD) CALL CALCGR
             IF(TYPMD.EQ.1) THEN
C           write md energies, etc. to .dat
              CALL MDRPREF(I)
             ELSE
              CALL MDRPRAB(I)
             ENDIF
C
             IF(TYPMD.EQ.1) THEN
              WRITE(IP,*) 'ROT.KE TEMP = ',TEMPTR(EKINR)
              WRITE(IP,*) 'TRANS.KE TEMP = ',TEMPTR(EKINT)
             ENDIF
                 JT = 0
                 ENDIF
            JT = JT + 1
C
C        *** every kevery timesteps: ***
C        ***     write configuration ***
C        ***     write g(r)s         ***
C
          IF ( KT .EQ. KEVERY ) THEN
C        write md configurations, velocities, grad, torq, com
         IF(TYPMD.EQ.1) THEN
              IF(MASWRK) CALL MDWREF(I,INTMD)
              IF(MASWRK) CALL TIMIT(1)
         ELSE
          CALL MDWRAB(I)
         ENDIF
             IF(PROD) CALL WRIGR
             KT = 0
          ENDIF
          KT = KT + 1
C
C
        ENDDO
C
        RETURN
C
C        *** error messages
C
 4999   FORMAT(/,1X,'***** PLEASE SPECIFY TYPE OF MD: ABINIT',
     *        ' OR EFP *****',/)
 5000         FORMAT(/,1X,'***** TIMESTEP (DT) TOO LARGE FOR STABLE',
     *               ' INTEGRATION *****',/)
 5050        FORMAT(/,1X,'***** ',I1,' ERROR(S) HAVE BEEN DETECTED',
     *               ' IN $MD. *****',/)
 5060        FORMAT(/,1X,'***** EXITING GAMESS. *****',/)
 6000        FORMAT(/,5X,'$MD OPTIONS',/,5X,11('-'),/,/,
     *         5X,'MDTYP  = ',A8,  9X,'MDINT  = ',A8,  /,
     *         5X,'BATHT  = ',F8.2,9X,'BOXL   = ',F8.4,/,
     *         5X,'DELR   = ',F8.2,9X,'DT     = ',E8.2,/,
     *         5X,'JEVERY = ',I8,  9X,'KEVERY = ',I8,  /,
     *         5X,'RETEMP = ',I8,  9X,'NSTEPS = ',I8,  /,
     *         5X,'NVT    = ',L8,  9X,'READ   = ',L8,  /,
     *         5X,'MBT    = ',L8,  9X,'MBR    = ',L8,  /,
     *         5X,'QRAND  = ',L8,  9X,'PROD   = ',L8,  /,
     *         5X,'EQUIL  = ',L8,  9X,'PBCOUT = ',L8     )
 6010 FORMAT(1X,'*** STARTING MD SIMULATION WITH ',I5,' FRAGMENTS ***'/)
 6011 FORMAT(1X,'*** STARTING MD SIMULATION WITH ',I5,' ATOMS ***'/)
 6020 FORMAT(5X,'INITIAL POTENTIAL ENERGY = ',F18.10,' AU')
 6030 FORMAT(5X,'INITIAL POTENTIAL ENERGY = ',F18.10,' KCAL/MOL',/)
      END
C
C-XXX-C*MODULE MDDRVR   *DECK TEMP1MD
C-XXX-      DOUBLE PRECISION FUNCTION TEMP1MD(EK)
C-XXX-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-XXX-      PARAMETER (MXFRG=1050, MXATM=2000,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-C
C-XXX-      DOUBLE PRECISION MDTYP
C-XXX-C
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
C-XXX-     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C-XXX-      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
C-XXX-C
C-XXX-      DATA EFP /8HEFP     /
C-XXX-C
C-XXX-C       tempmd = 2.0d+00 * ek / (6.0d+00*nfrg-3.0d+00)
C-XXX-C     *                 / 1.987216d-03
C-XXX-C--pre-IBM
C-XXX-C          temp1md = 2.0d+00 * ek / (6.0d+00*nfrg-6.0d+00)
C-XXX-C     *                 / 1.987216d-03
C-XXX-C
C-XXX-      IF(MDTYP.EQ.EFP) THEN
C-XXX-         TEMP1MD = 2.0D+00 * EK / (6.0D+00*NFRG) / 1.987216D-03
C-XXX-      ELSE
C-XXX-         TEMP1MD = 2.0D+00 * EK / (6.0D+00*NAT)  / 1.987216D-03
C-XXX-      ENDIF
C-XXX-      RETURN
C-XXX-      END
C
C*MODULE MDDRVR   *DECK TEMPTR
      DOUBLE PRECISION FUNCTION TEMPTR(EK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050,
     *             MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
C       tempmd = 2.0d+00 * ek / (6.0d+00*nfrg-3.0d+00)
C     *                 / 1.987216d-03
C--pre-IBM
C        temptr = 2.0d+00 * ek / (3.0d+00*nfrg-3.0d+00)
C     *                 / 1.987216d-03
C
      TEMPTR = 2.0D+00 * EK / (3.0D+00*NFRG) / 1.987216D-03
      RETURN
      END
C
C*MODULE MDMAIN   *DECK ENEDAB
      SUBROUTINE ENEDAB
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /FUNCT / E, EG(3*MXATM)
C
C        ***************************************
C        *** Calculate energy and forces     ***
C        ***************************************
C
      WRITE(6,*) 'IN ENEDAB...'
      CALL GRADX
      F2F = -AU2KCAL / EK2KCAL
      DO I = 1, NAT*3
         EG(I) = EG(I)*F2F
      ENDDO
      RETURN
      END
C*MODULE MDMAIN   *DECK ENEDEF
      SUBROUTINE ENEDEF
C       *********************************************************
C       *** Calculate energy and forces                       ***
C       *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFGPT=12000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /EFPOTD/ METHOD,INABIO,MOVE
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
      IF (METHOD.EQ.6) THEN
         CALL EFCM
         CALL SPCEG
         CALL SPCTGT
      ELSE
         CALL EFSP
         CALL EFGRAD
      ENDIF
C
      F2F = -AU2KCAL / EK2KCAL
C
C TESTING GRADIENT AND TORQUE
C        do j=1,nfrg
C           write(6,9035) j,deft(1,j),deft(2,j),deft(3,j)
C        enddo
C 9035  FORMAT(2X,'fragment ',2X,i6,2X,'GRADIENT  ',4X,3F18.9)
C
C        do j=1,nfrg
C           write(6,9036) j,torq(1,j),torq(2,j),torq(3,j)
C        enddo
C 9036  FORMAT(2X,'fragment ',2X,i6,2X,'TORQUE  ',4X,3F18.9)
C
       DO I = 1, NFRG
          DEFT(1,I) = F2F * DEFT(1,I)
          DEFT(2,I) = F2F * DEFT(2,I)
          DEFT(3,I) = F2F * DEFT(3,I)
          TORQ(1,I) = F2F * TORQ(1,I)
          TORQ(2,I) = F2F * TORQ(2,I)
          TORQ(3,I) = F2F * TORQ(3,I)
       ENDDO
       RETURN
       END
C
C*MODULE MDMAIN   *DECK MDMOVE(INTMD)
      SUBROUTINE MDMOVE(INTMD)
C        *********************************************************
C        *** Integrate equations of motion using a leap-frog   ***
C        *** algorithm and quaternion                          ***
C        *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXPT=2000, MXFGPT=12000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MXAO=8192, MXSH=5000, MXGTOT=20000)
      DOUBLE PRECISION MDTYP
      INTEGER INTMD
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /FMCOM / XX(1)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C
C        write(6,*) 'vx(1) = ',vx(1)
C        call mdwrit(1)
        EKINT = 0.0D+00
        EKINR = 0.0D+00
        NMST = 0
        NPST = 0
        NDPST = 0
        NRST = 0
        NATST = 0
        NORST = 0
        NAOST = 0
        NDST = 0
C
C
      MXBF=0
      DO I=1,NFRG
         MXBF=MAX(MXBF,NPBF(I))
      END DO
      CALL VALFM(LOADFM)
      LPROVEC = LOADFM  + 1
      LCANVEC = LPROVEC + MXBF*NTMO
      LAST    = LCANVEC + MXBF*NTAO
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
      LENPV=MXBF*NTMO
      IF (NFRG.GT.0.AND.NTMO.GT.0)
     *   CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
      LENCV=MXBF*NTAO
      IF (NFRG.GT.0.AND.NTAO.GT.0)
     *   CALL DAREAD(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
C
        DT2 = 0.5D+00*DT
        DEL = 0.0D+00
        DO I = 1, NFRG
C
C          write(6,*) 'torq(1,i),(2,i),(3,i) = ',torq(1,i),
C     *                torq(2,i),torq(3,i)
C   Move angular momentum to time t-->
C             oxi(t)= ox(i)(t-1/2dt) + 1/2dt*torq(t)
C
          OXI = OX(I) + DT2 * TORQ(1,I)
          OYI = OY(I) + DT2 * TORQ(2,I)
          OZI = OZ(I) + DT2 * TORQ(3,I)
C
C   Get rotation matrix at time t--> A(t) from quat(i)(t)
          AXX = QW(I) ** 2 + QX(I) ** 2 - QY(I) ** 2 - QZ(I) ** 2
          AXY = 2.0D+00 * ( QX(I) * QY(I) + QW(I) * QZ(I) )
          AXZ = 2.0D+00 * ( QX(I) * QZ(I) - QW(I) * QY(I) )
          AYX = 2.0D+00 * ( QX(I) * QY(I) - QW(I) * QZ(I) )
          AYY = QW(I) ** 2 - QX(I) ** 2 + QY(I) ** 2 - QZ(I) ** 2
          AYZ = 2.0D+00 * ( QY(I) * QZ(I) + QW(I) * QX(I) )
          AZX = 2.0D+00 * ( QX(I) * QZ(I) + QW(I) * QY(I) )
          AZY = 2.0D+00 * ( QY(I) * QZ(I) - QW(I) * QX(I) )
          AZZ = QW(I) ** 2 - QX(I) ** 2 - QY(I) ** 2 + QZ(I) ** 2
C
C   Save quaternions at time t
          QWSAV = QW(I)
          QXSAV = QX(I)
          QYSAV = QY(I)
          QZSAV = QZ(I)
C
C   Convert angular mom to body-fixed form--to angular vel-->
C   omxi(t) = A(t) * ang mom oxi(t)
          OMXI = ( AXX * OXI + AXY * OYI + AXZ * OZI ) / AMOM(1,ISET(I))
          OMYI = ( AYX * OXI + AYY * OYI + AYZ * OZI ) / AMOM(2,ISET(I))
          OMZI = ( AZX * OXI + AZY * OYI + AZZ * OZI ) / AMOM(3,ISET(I))
C
C          write(6,*) 'i, ang vel = ',omxi*1.0d-12,omyi*1.0d-12,
C     *                               omzi*1.0d-12
C   Kinetic energy at time t
      EKINR = EKINR + AMOM(1,ISET(I)) * OMXI**2 +
     *        AMOM(2,ISET(I)) * OMYI**2 + AMOM(3,ISET(I)) * OMZI**2
C
C   Obtain time-deriv of quat. and advance to time t+1/2dt-->
C   qwli-q'(t) = 1/2*q(t)*ang. vel omxi(t)
C   Guess at q(t+1/2dt)=q(t) + 1/2dt*q'(t)
          QW1I = ( - QX(I)*OMXI - QY(I)*OMYI - QZ(I)*OMZI )*0.5D+00
          QX1I = (   QW(I)*OMXI - QZ(I)*OMYI + QY(I)*OMZI )*0.5D+00
          QY1I = (   QZ(I)*OMXI + QW(I)*OMYI - QX(I)*OMZI )*0.5D+00
          QZ1I = ( - QY(I)*OMXI + QX(I)*OMYI + QW(I)*OMZI )*0.5D+00
C
C          write(6,*) '1st time deriv; ',qw1i,qx1i,qy1i,qz1i
C--derivative constraint?
          DEL = QW(I)*QW1I+QX(I)*QX1I+QY(I)*QY1I+QZ(I)*QZ1I
          QW1I = QW1I - DEL*QW(I)
          QX1I = QX1I - DEL*QX(I)
          QY1I = QY1I - DEL*QY(I)
          QZ1I = QZ1I - DEL*QZ(I)
C
          QWI = QW(I) + DT2 * QW1I
          QXI = QX(I) + DT2 * QX1I
          QYI = QY(I) + DT2 * QY1I
          QZI = QZ(I) + DT2 * QZ1I
C
C        the following gives a little (just a little) more stability
C
C   Renormalize new quat. at t+1/2dt
          QLI = 1.0D+00/SQRT(QWI**2 + QXI**2 + QYI**2 + QZI**2)
          QWI = QWI * QLI
          QXI = QXI * QLI
          QYI = QYI * QLI
          QZI = QZI * QLI
C
C   Form new rotation matrix at t+1/2dt
          AXX = QWI ** 2 + QXI ** 2 - QYI ** 2 - QZI ** 2
          AXY = 2.0D+00 * ( QXI * QYI + QWI * QZI )
          AXZ = 2.0D+00 * ( QXI * QZI - QWI * QYI )
          AYX = 2.0D+00 * ( QXI * QYI - QWI * QZI )
          AYY = QWI ** 2 - QXI ** 2 + QYI ** 2 - QZI ** 2
          AYZ = 2.0D+00 * ( QYI * QZI + QWI * QXI )
          AZX = 2.0D+00 * ( QXI * QZI + QWI * QYI )
          AZY = 2.0D+00 * ( QYI * QZI - QWI * QXI )
          AZZ = QWI ** 2 - QXI ** 2 - QYI ** 2 + QZI ** 2
C
C   Move ang. mom from t-1/2dt to t+1/2dt and store-->
C   ox(i)(t+1/2dt) = ox(i)(t-1/2dt) + dt*torq(t)
          OX(I) = OX(I) + DT*TORQ(1,I)
          OY(I) = OY(I) + DT*TORQ(2,I)
          OZ(I) = OZ(I) + DT*TORQ(3,I)
C
C   Convert ang. mom to body-fixed form to get ang. vel at time
C   t+1/2dt-->    omxi(t+1/2dt) = A(t+1/2dt)*ox(i)(t+1/2dt)
C
       OMXI=(AXX*OX(I) + AXY*OY(I) + AXZ*OZ(I))/AMOM(1,ISET(I))
       OMYI=(AYX*OX(I) + AYY*OY(I) + AYZ*OZ(I))/AMOM(2,ISET(I))
       OMZI=(AZX*OX(I) + AZY*OY(I) + AZZ*OZ(I))/AMOM(3,ISET(I))
C
C          write(6,*) 'i, ang vel = ',omxi*1.0d-12,omyi*1.0d-12,
C     *                               omzi*1.0d-12
C   Obtain time-deriv of quat. at t+1/2dt and advance to time t+dt-->
C   qwli-q'(t+1/2dt) = 1/2*q(t+1/2dt)*ang. vel omxi(t+1/2dt)
C   Guess at q(t+dt)=q(t) + dt*q'(t+1/2dt)
C
       QW1I=(-QXI*OMXI - QYI*OMYI - QZI*OMZI )*0.5D+00
       QX1I=( QWI*OMXI - QZI*OMYI + QYI*OMZI )*0.5D+00
       QY1I=( QZI*OMXI + QWI*OMYI - QXI*OMZI )*0.5D+00
       QZ1I=(-QYI*OMXI + QXI*OMYI + QWI*OMZI )*0.5D+00
C
          WRITE(6,*) 'I,TIME DERIV = ',QW1I,QX1I,QY1I,QZ1I
C
          QW(I) = QW(I) + DT*QW1I
          QX(I) = QX(I) + DT*QX1I
          QY(I) = QY(I) + DT*QY1I
          QZ(I) = QZ(I) + DT*QZ1I
C
C   Renormalize new quat. at t+dt
          QLI = 1.0D+00/SQRT(QW(I)**2 + QX(I)**2 + QY(I)**2 + QZ(I)**2)
C          write(6,*) 'renorm factor ',qli
C
          QW(I) = QW(I) * QLI
          QX(I) = QX(I) * QLI
          QY(I) = QY(I) * QLI
          QZ(I) = QZ(I) * QLI
C
C          write(6,*) 'qw,x,y,z = ',qw(i),qx(i),qy(i),qz(i)
C   Save vel. at time t-1/2dt
          VXI = VX(I)
          VYI = VY(I)
          VZI = VZ(I)
C
C   Move linear vel. from t-1/2dt to t+1/2dt and save-->
C   vx(i)(t+1/2dt) = vxi(t-1/2dt) + dt*force(t)/mass
          VX(I) = VXI + DT*DEFT(1,I) * TOTMAI(ISET(I))
          VY(I) = VYI + DT*DEFT(2,I) * TOTMAI(ISET(I))
          VZ(I) = VZI + DT*DEFT(3,I) * TOTMAI(ISET(I))
C
C   Calc. current velocity (at t) to evaluate energy, etc.-->
C   vxi(t) = 1/2*(vxi(t-1/2dt)+vx(i)(t+1/2dt))
          VXI = 0.5D+00 * ( VXI + VX(I) )
          VYI = 0.5D+00 * ( VYI + VY(I) )
          VZI = 0.5D+00 * ( VZI + VZ(I) )
C
C   Calc. KEt at time t
      EKINT = EKINT + TOTMAS(ISET(I))*(VXI**2 + VYI**2 + VZI**2)
C
C   Save new vel. at t+1/2dt; note that new positions at t+dt-->
C   r(t+dt) = r(t) + dt*v(t+1/2dt)
          CENTDX = DT*VX(I)
          CENTDY = DT*VY(I)
          CENTDZ = DT*VZ(I)
C
C          write(6,*) 'centdx,y,z = ',centdx,centdy,centdz
C  Do rotation and translation to get ecp to t+dt
C  Note: we need to call efcm before coming back
C        here to get newest COM coords
C
          CALL MDDSPQ(I,CENTDX,CENTDY,CENTDZ,QWSAV,QXSAV,QYSAV,QZSAV,
     *                QW(I),QX(I),QY(I),QZ(I),
     *                NMST,NPST,NDPST,NRST,NDST,NATST,NORST,
     *                XX(LPROVEC),NAOST,XX(LCANVEC),MXBF,INTMD)
C
          NMST  = NMST  + NMPTS(I)
          NPST  = NPST  + NPPTS(I)
          NDPST = NDPST + NDPPTS(I)
          NRST  = NRST  + NRPTS(I)
          NDST  = NDST  + NDPTS(I)
          NATST = NATST + NATEF(ISET(I))
          NORST = NORST + NORB(I)
          NAOST = NAOST + NAO(ISET(I))
C
        ENDDO
C
      IF (NFRG.GT.0.AND.NTMO.GT.0)
     *   CALL DAWRIT(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
      IF (NFRG.GT.0.AND.NTAO.GT.0)
     *   CALL DAWRIT(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
C
      CALL RETFM(NEED)
C
        EKINT = 0.5D+00 * EK2KCAL * EKINT
        EKINR = 0.5D+00 * EK2KCAL * EKINR
      CALL RETFM(NEED)
        RETURN
        END
C*MODULE MDMAIN   *DECK MDMOVA(EG,NAT,C)
      SUBROUTINE MDMOVA(EG,NAT,C)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050)
C
      DOUBLE PRECISION MDTYP
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /HMASS/  TOTALM,ZMASSI(MXATM)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
C
      DIMENSION EG(3,NAT),C(3,NAT)
C
C       *********************************************************
C       *** Integrate equations of motion using a leap-frog   ***
C       *** algorithm and quaternions                         ***
C       *********************************************************
C
      EKINT = 0.0D+00
      TOTALM = 0.0D+00
      DO I = 1,NAT
          ZMASSI(I) = 1.0D+00 / ZMASS(I)
C
          VXI = VX(I)
          VYI = VY(I)
          VZI = VZ(I)
C
          VX(I) = VXI + EG(1,I) * ZMASSI(I)
          VY(I) = VYI + EG(2,I) * ZMASSI(I)
          VZ(I) = VZI + EG(3,I) * ZMASSI(I)
C
          VXI = 0.5D+00 * ( VXI + VX(I) )
          VYI = 0.5D+00 * ( VYI + VY(I) )
          VZI = 0.5D+00 * ( VZI + VZ(I) )
C
          EKINT = EKINT + VXI**2 + VYI**2 + VZI**2
C
          C(1,I) = C(1,I) + VX(I)
          C(2,I) = C(2,I) + VY(I)
          C(3,I) = C(3,I) + VZ(I)
C
          TOTALM = TOTALM + ZMASS(I)
        ENDDO
         EKINT = 0.5D+00 * EK2KCAL * TOTALM * EKINT
C
        RETURN
        END
C
C*MODULE MDINIT   *DECK MDCONV
      SUBROUTINE MDCONV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
C
C        *********************************************************
C        *** setup convenient conversion factors for MD        ***
C        *********************************************************
C        conversion factors taken (where possible) from NIST 1986 values
C
      AU2KCAL = 4.359816D-18 / 1.380658D-23 * 1.98721984D-03
      AU2ANG = 0.52917724924D+00
      EK2KCAL = 2.39005738055577D-27 * AU2ANG * AU2ANG
      AU2SEC = 2.41889D-17
C--hmn (au--kg--amu)
      AU2AMU = 9.1091D-31 / 1.6605402D-27
      RETURN
      END
C
C*MODULE MDINIT   *DECK MODSET
      SUBROUTINE MODSET
C        *******************************************
C        *** setup water model parameters for MD ***
C        *******************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MDTYP
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      PARAMETER (MXPT=2000, MXFGPT=12000,MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (AVOGAD=6.022045D+23,ZERO=0.0D+00)
C
      CHARACTER*8 FRGNME
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        body-fixed moments of inertia are derived from fragment #1
C        body-fixed frame is as follows:
C               z-axis is C2 axis (pointing away from hydrogens)
C               x-axis is pointing out of plane of hydrogens
C               y-axis points in plane towards H2
C        (forming a right-handed coordinate system)
C
        ROH1X = EFC(1,2) - EFC(1,1)
        ROH1Y = EFC(2,2) - EFC(2,1)
        ROH1Z = EFC(3,2) - EFC(3,1)
        ROH2X = EFC(1,3) - EFC(1,1)
        ROH2Y = EFC(2,3) - EFC(2,1)
        ROH2Z = EFC(3,3) - EFC(3,1)
        ZX = -(ROH1X + ROH2X)
        ZY = -(ROH1Y + ROH2Y)
        ZZ = -(ROH1Z + ROH2Z)
        ZL = SQRT(ZX*ZX+ZY*ZY+ZZ*ZZ)
        ZX = ZX / ZL
        ZY = ZY / ZL
        ZZ = ZZ / ZL
        XX = ZY * ROH1Z - ZZ * ROH1Y
        XY = ZZ * ROH1X - ZX * ROH1Z
        XZ = ZX * ROH1Y - ZY * ROH1X
        XL = SQRT(XX*XX+XY*XY+XZ*XZ)
        XX = XX / XL
        XY = XY / XL
        XZ = XZ / XL
        YX = ZY * XZ - ZZ * XY
        YY = ZZ * XX - ZX * XZ
        YZ = ZX * XY - ZY * XX
        YL = SQRT(YX*YX+YY*YY+YZ*YZ)
        YX = YX / YL
        YY = YY / YL
        YZ = YZ / YL
C
           L=0
C
       FMASAL=ZERO
      DO 5 K=1,NFRG
        TOTMAS(ISET(K)) = FRGMAS(K)
        TOTMAI(ISET(K)) = 1.0D+00 / TOTMAS(ISET(K))
        AMOM(1,ISET(K)) = 0.0D+00
        AMOM(2,ISET(K)) = 0.0D+00
        AMOM(3,ISET(K)) = 0.0D+00
C
       DO I = 1, NMPTS(K)
         DRX(I+L) = (EFC(1,I) - EFCENT(1,1))*XX
     *          + (EFC(2,I) - EFCENT(2,1))*XY
     *          + (EFC(3,I) - EFCENT(3,1))*XZ
          DRY(I+L) = (EFC(1,I) - EFCENT(1,1))*YX
     *          + (EFC(2,I) - EFCENT(2,1))*YY
     *          + (EFC(3,I) - EFCENT(3,1))*YZ
          DRZ(I+L) = (EFC(1,I) - EFCENT(1,1))*ZX
     *          + (EFC(2,I) - EFCENT(2,1))*ZY
     *          + (EFC(3,I) - EFCENT(3,1))*ZZ
C
          AMOM(1,ISET(K)) = AMOM(1,ISET(K)) + FMASS(I+L) *
     *                      (DRY(I+L)**2 + DRZ(I+L)**2)
          AMOM(2,ISET(K)) = AMOM(2,ISET(K)) + FMASS(I+L) *
     *                      (DRX(I+L)**2 + DRZ(I+L)**2)
          AMOM(3,ISET(K)) = AMOM(3,ISET(K)) + FMASS(I+L) *
     *                     (DRX(I+L)**2 + DRY(I+L)**2)
C
C
        ENDDO
C
         FMASAL=FMASAL+FRGMAS(K)
C
          L = L + NMPTS(K)
    5       CONTINUE
C
C
         WMAS=18.01056D+00
        BOXLAN = BOXL*AU2ANG
          IF(NDFRG.EQ.1 .AND. FRGMAS(1).EQ.WMAS) THEN
C       Note that watmas = mass of water molecule in g/mol, simden =
C       density of simulation system assuming all molecules are like
C       molecule #1. Must change volume from bohr^3 to Ang^3 to cm^3.
        WATMAS = ZERO
        SIMDEN = ZERO
        SIMDEN1 = ZERO
        DO I = 1,3
         WATMAS = WATMAS + FMASS(I)
        ENDDO
C
C
        SIMDEN=(NFRG/AVOGAD*WATMAS)/(BOXLAN*BOXLAN*BOXLAN
     *          *1.0D-8*1.0D-8*1.0D-8)
        IF(MASWRK) WRITE(IW,2000) NFRG,BOXL*AU2ANG,SIMDEN
C
        VOL1 = ZERO
        BOXL1 = ZERO
        VOL1 = (NFRG/AVOGAD*WATMAS)/(1.0D-8)**3/0.99697D+00
        CUBERT=1.0D+00/3.0D+00
        BOXL1 = VOL1**CUBERT
        IF(MASWRK) WRITE(IW,2010) NFRG,BOXL1
C
        ELSE
        SIMDEN1=(FMASAL/AVOGAD)/(BOXLAN*BOXLAN*BOXLAN
     *          *1.0D-8*1.0D-8*1.0D-8)
              IF(MASWRK) WRITE(IW,2001) NFRG,BOXL*AU2ANG,SIMDEN1
        END IF
C
 2000   FORMAT(1X,'***MD SIMULATION WITH ',I5,
     *         ' WATER MOLECULES',/,1X,'***BOXLENGTH = ',F8.4,
     *         ' ANGSTROMS',/,1X,'***NVT DENSITY = ',F8.4,' G/ML')
 2001   FORMAT(1X,'***MD SIMULATION WITH ',I5,
     *         ' MOLECULES',/,1X,'***BOXLENGTH = ',F8.4,
     *         ' ANGSTROMS',/,1X,'***NVT DENSITY = ',F8.4,' G/ML')
 2010   FORMAT(1X,'***FOR DENSITY OF 0.997 G/ML (300K EXPT) WITH ',I5,
     *         ' WATER MOLECULES, THEORETICAL BOXLENGTH = ',F8.4,
     *         ' ANGSTROMS',/)
C
        RETURN
        END
C
C*MODULE MDINIT   *DECK CLRTVE
      SUBROUTINE CLRTVE
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
C        ****************************
C        *** setup arrays for MD  ***
C        ****************************
C
      DO I = 1, NFRG
         VX(I) = ZERO
         VY(I) = ZERO
         VZ(I) = ZERO
      ENDDO
      RETURN
      END
C
C*MODULE MDINIT   *DECK CLRMOM
      SUBROUTINE CLRMOM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
C       *********************************************************
C       *** setup rotational arrays for MD                    ***
C       *********************************************************
C
      DO I = 1, NFRG
         OX(I) = ZERO
         OY(I) = ZERO
         OZ(I) = ZERO
      ENDDO
      RETURN
      END
C
C*MODULE MDINIT   *DECK CLRVEL
      SUBROUTINE CLRVEL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
C       *********************************************************
C       *** setup rotational arrays for MD                    ***
C       *********************************************************
C
      DO I = 1, NFRG
         WX(I) = ZERO
         WY(I) = ZERO
         WZ(I) = ZERO
      ENDDO
      RETURN
      END
C
C-XXX-c*module mdinit   *deck cart2q
C-XXX-      subroutine cart2q
C-XXX-c        *********************************************************
C-XXX-c        *** convert cartesian coordinates to quaternions      ***
C-XXX-c        *** (in an ironic twist of programming savviness,     ***
C-XXX-c        ***  this is a fortran port of my original perl       ***
C-XXX-c        ***  script)         ***
C-XXX-c        *********************************************************
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-c
C-XXX-      PARAMETER (MXFGPT=12000, MXFRG=1050,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-c
C-XXX-      CHARACTER*8 FRGNME
C-XXX-c
C-XXX-      logical        found
C-XXX-c
C-XXX-      dimension d1x(3), d1y(3), d1z(3)
C-XXX-c
C-XXX-      common /atmdat/ vx(mxfrg), vy(mxfrg), vz(mxfrg),
C-XXX-     *                qw(mxfrg), qx(mxfrg), qy(mxfrg), qz(mxfrg),
C-XXX-     *                ox(mxfrg), oy(mxfrg), oz(mxfrg), qw1(mxfrg),
C-XXX-     *                qx1(mxfrg),qy1(mxfrg),qz1(mxfrg),qw2(mxfrg),
C-XXX-     *                qx2(mxfrg),qy2(mxfrg),qz2(mxfrg)
C-XXX-      common /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),
C-XXX-     *                EFATRM(MXFGPT),
C-XXX-     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),
C-XXX-     *                EFBTRM2(MXFGPT),
C-XXX-     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
C-XXX-     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
C-XXX-      common /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
C-XXX-     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),
C-XXX-     *                FRGMI(6,MXFRG),ATORQ(3,MXFRG)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      common /iofile/ ir,iw,ip,is,ipk,idaf,nav,ioda(950)
C-XXX-c
C-XXX-      data toler /1.0d-4/
C-XXX-c
C-XXX-      k = 0
C-XXX-      do i = 1, nfrg
C-XXX-c
C-XXX-          do j = 1,3
C-XXX-            d1x(j) = efc(1,j+k) - efcent(1,i)
C-XXX-            d1y(j) = efc(2,j+k) - efcent(2,i)
C-XXX-            d1z(j) = efc(3,j+k) - efcent(3,i)
C-XXX-          enddo
C-XXX-c
C-XXX-          wzx = -(d1x(2) + d1x(3))
C-XXX-          wzy = -(d1y(2) + d1y(3))
C-XXX-          wzz = -(d1z(2) + d1z(3))
C-XXX-          wl = sqrt(wzx*wzx+wzy*wzy+wzz*wzz)
C-XXX-          wzx = wzx / wl
C-XXX-          wzy = wzy / wl
C-XXX-          wzz = wzz / wl
C-XXX-c
C-XXX-          wxx = wzy * d1z(2) - wzz * d1y(2)
C-XXX-          wxy = wzz * d1x(2) - wzx * d1z(2)
C-XXX-          wxz = wzx * d1y(2) - wzy * d1x(2)
C-XXX-          wl = sqrt(wxx*wxx+wxy*wxy+wxz*wxz)
C-XXX-          wxx = wxx / wl
C-XXX-          wxy = wxy / wl
C-XXX-          wxz = wxz / wl
C-XXX-c
C-XXX-          wyx = wzy * wxz - wzz * wxy
C-XXX-          wyy = wzz * wxx - wzx * wxz
C-XXX-          wyz = wzx * wxy - wzy * wxx
C-XXX-          wl = sqrt(wyx*wyx+wyy*wyy+wyz*wyz)
C-XXX-          wyx = wyx / wl
C-XXX-          wyy = wyy / wl
C-XXX-          wyz = wyz / wl
C-XXX-c
C-XXX-          rat = (wxz + wzx) / (wxy + wyx)
C-XXX-          wqz = sqrt(rat * ( wyz + wzy ) * 0.25d+00)
C-XXX-          wqy = wqz / rat
C-XXX-          wqx = 0.25d+00 * ( wxz + wzx ) / wqz
C-XXX-          wqw = sqrt(1.0d+00 - wqx*wqx - wqy*wqy - wqz*wqz)
C-XXX-c
C-XXX-          found = .false.
C-XXX-c
C-XXX-          do l = 1, 4
C-XXX-            if ( l .eq. 1 ) then
C-XXX-              wqwl = wqw
C-XXX-              wqzl = wqz
C-XXX-            else if ( l .eq. 2 ) then
C-XXX-              wqwl = wqw
C-XXX-              wqzl = -wqz
C-XXX-            else if ( l .eq. 3 ) then
C-XXX-              wqwl = -wqw
C-XXX-              wqzl = wqz
C-XXX-            else
C-XXX-              wqwl = -wqw
C-XXX-              wqzl = -wqz
C-XXX-            endif
C-XXX-c
C-XXX-            wqyl = wqzl / rat
C-XXX-            wqxl = 0.25d+00 * ( wxz + wzx ) / wqzl
C-XXX-c
C-XXX-            axx = wqwl ** 2 + wqxl ** 2 - wqyl ** 2 - wqzl ** 2
C-XXX-            axy = 2.0d+00 * ( wqxl * wqyl + wqwl * wqzl )
C-XXX-            axz = 2.0d+00 * ( wqxl * wqzl - wqwl * wqyl )
C-XXX-            ayx = 2.0d+00 * ( wqxl * wqyl - wqwl * wqzl )
C-XXX-            ayy = wqwl ** 2 - wqxl ** 2 + wqyl ** 2 - wqzl ** 2
C-XXX-            ayz = 2.0d+00 * ( wqyl * wqzl + wqwl * wqxl )
C-XXX-            azx = 2.0d+00 * ( wqxl * wqzl + wqwl * wqyl )
C-XXX-            azy = 2.0d+00 * ( wqyl * wqzl - wqwl * wqxl )
C-XXX-            azz = wqwl ** 2 - wqxl ** 2 - wqyl ** 2 + wqzl ** 2
C-XXX-c
C-XXX-            error = (axx-wxx)**2+(axy-wxy)**2+(axz-wxz)**2
C-XXX-     *            + (ayx-wyx)**2+(ayy-wyy)**2+(ayz-wyz)**2
C-XXX-     *            + (azx-wzx)**2+(azy-wzy)**2+(azz-wzz)**2
C-XXX-c
C-XXX-            if ( error .lt. toler ) then
C-XXX-              found = .true.
C-XXX-              goto 100
C-XXX-            endif
C-XXX-          enddo
C-XXX-c
C-XXX- 100          if ( found ) then
C-XXX-            qw(i) = wqwl
C-XXX-            qx(i) = wqxl
C-XXX-            qy(i) = wqyl
C-XXX-            qz(i) = wqzl
C-XXX-          else
C-XXX-            write(iw,1000) i
C-XXX-            call abrt
C-XXX-          endif
C-XXX-          k = k + nmpts(i)
C-XXX-       enddo
C-XXX-c
C-XXX-       return
C-XXX- 1000  format(/,1x,'*** ERROR CONVERTING CART. COORDINATES TO ',
C-XXX-     *             'QUATERNIONS: FRAGMENT ',i8)
C-XXX-       end
C
C*MODULE MDINIT   *DECK INITGR
      SUBROUTINE INITGR
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (NBIN=1000, ZERO=0.0D+00)
C
      DOUBLE PRECISION MDTYP
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /GOFR  / GOFR(NBIN,3), DELR, NUMG, NUMSAM
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
C
C        *********************************************************
C        *** initialize g(r) arrays, g(r) parameters      ***
C        *********************************************************
C
      NUMG = INT( BOXL * AU2ANG * 0.5D+00 / DELR ) + 1
      IF ( NUMG .GT. NBIN ) THEN
        NUMG = NBIN
      ENDIF
      DO I = 1, 3
        DO J = 1, NUMG
          GOFR(J, I) = ZERO
        ENDDO
      ENDDO
      NUMSAM = 0
      RETURN
      END
C
C*MODULE MDINIT   *DECK RANDQ
      SUBROUTINE RANDQ
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL DEBUG
C
      DIMENSION ITIMST(6)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
C
      DEBUG=.TRUE.
      IF (DEBUG) THEN
         IX = 666000666
      ELSE
      CALL TMDATE(ITIMST)
      IX=0
      LENTIM=6
      IF(NWDVAR.EQ.1) LENTIM=3
      DO I=1,LENTIM
        IX=IX+ITIMST(I)
      ENDDO
      IX = IABS(IX)
      ENDIF
      IDUM = -IX
C
      DO I = 1, NFRG
C
         CALL RAND3(IDUM, RANM1)
         CALL RAND3(IDUM, RANM2)
         CALL RAND3(IDUM, RANM3)
         CALL RAND3(IDUM, RANM4)
          WX = 2.0D+00*RANM1-1.0D+00
          WY = 2.0D+00*RANM2-1.0D+00
          WZ = 2.0D+00*RANM3-1.0D+00
          WW = 2.0D+00*RANM4-1.0D+00
          RAT = SQRT(WX*WX + WY*WY + WZ*WZ + WW*WW)
          WXN = WX / RAT
          WYN = WY / RAT
          WZN = WZ / RAT
          WWN = WW / RAT
C
          QW(I) = WWN
          QX(I) = WXN
          QY(I) = WYN
          QZ(I) = WZN
      ENDDO
      RETURN
      END
C
C-XXX-c*module mdinit   *deck comfy
C-XXX-      subroutine comfy
C-XXX-c
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-      double precision max2d
C-XXX-C
C-XXX-      PARAMETER (MXFGPT=12000,MXFRG=1050,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-c
C-XXX-      common /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
C-XXX-     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),
C-XXX-     *                FRGMI(6,MXFRG),ATORQ(3,MXFRG)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      common /SYSCET/ SYSMAS,SYSCNT(3)
C-XXX-C
C-XXX-      DATA ZERO/0.0D+00/
C-XXX-c
C-XXX-      max2d = zero
C-XXX-C      maxd= zero
C-XXX-C      comfyl = zero
C-XXX-c
C-XXX-      do i=1,nfrg
C-XXX-        xd = syscnt(1) - efcent(1,i)
C-XXX-        yd = syscnt(2) - efcent(2,i)
C-XXX-        zd = syscnt(3) - efcent(3,i)
C-XXX-        r2d = xd*xd + yd*yd + zd*zd
C-XXX-c
C-XXX-        if(r2d .gt. max2d) then
C-XXX-          max2d = r2d
C-XXX-        endif
C-XXX-      enddo
C-XXX-c
C-XXX-C      maxd = sqrt(max2d)
C-XXX-C      comfyl = 2*maxd + 2.0d+00/au2ang
C-XXX-c
C-XXX-      return
C-XXX-      end
C
C-XXX-c*module mdinit   *deck gaussv
C-XXX-      SUBROUTINE GAUSSV(temp)
C-XXX-c
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-C--      double precision m2bohr
C-XXX-C
C-XXX-      PARAMETER ( MXFGPT=12000,MXFRG=1050,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-c
C-XXX-      common /atmdat/ vx(mxfrg), vy(mxfrg), vz(mxfrg),
C-XXX-     *                qw(mxfrg), qx(mxfrg), qy(mxfrg), qz(mxfrg),
C-XXX-     *                ox(mxfrg), oy(mxfrg), oz(mxfrg), qw1(mxfrg),
C-XXX-     *                qx1(mxfrg),qy1(mxfrg),qz1(mxfrg),qw2(mxfrg),
C-XXX-     *                qx2(mxfrg),qy2(mxfrg),qz2(mxfrg)
C-XXX-      common /convmd/ au2kcal, ek2kcal, au2ang, au2sec, au2amu
C-XXX-      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
C-XXX-     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),
C-XXX-     *                FRGMI(6,MXFRG),ATORQ(3,MXFRG)
C-XXX-      common /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-c
C-XXX-      DATA HALF,ONE,TWO/0.5d+00,1.0d+00,2.0d+00/
C-XXX-C
C-XXX-c         initialize RNG
C-XXX-c
C-XXX-      iseed = -43891265
C-XXX-      call ranth(iseed,rannum)
C-XXX-      iseed = 0
C-XXX-c
C-XXX-c      m2bohr=1.0d+00/5.2917706d-11
C-XXX-c      boltzk = 1.3806503d-23/1.66053886d-27*m2bohr*m2bohr
C-XXX-      boltzk = 0.83112120d+00/1.0d-12/1.0d-12/au2ang/au2ang
C-XXX-      TWOPI=TWO*ACOS(-ONE)
C-XXX-      DO 2 I=1,nfrg
C-XXX-          SIGMA=SQRT((BOLTZK*TEMP)/frgmas(I))
C-XXX-          CALL RANth(ISEED,R1)
C-XXX-          CALL RANth(ISEED,R2)
C-XXX-          FACT=SQRT(-TWO*LOG(R1))
C-XXX-          RANN1=FACT*COS(TWOPI*R2)
C-XXX-          RANN2=FACT*SIN(TWOPI*R2)
C-XXX-          CALL RANth(ISEED,R1)
C-XXX-          CALL RANth(ISEED,R2)
C-XXX-          FACT=SQRT(-TWO*LOG(R1))
C-XXX-          RANN3=FACT*COS(TWOPI*R2)
C-XXX-          CALL RANth(ISEED,R1)
C-XXX-          CALL RANth(ISEED,R2)
C-XXX-          CALL RANth(ISEED,R3)
C-XXX-          IF (R1.GT.HALF)RANN1=-RANN1
C-XXX-          IF (R2.GT.HALF)RANN2=-RANN2
C-XXX-          IF (R3.GT.HALF)RANN3=-RANN3
C-XXX-          VX(I)=SIGMA*RANN1
C-XXX-          VY(I)=SIGMA*RANN2
C-XXX-          VZ(I)=SIGMA*RANN3
C-XXX-C    1   CONTINUE
C-XXX-    2 CONTINUE
C-XXX-c       Now v units are bohr/s
C-XXX-      RETURN
C-XXX-      END
C
C-XXX-c*module mdinit   *deck ranth
C-XXX-      SUBROUTINE RANth(IDUM,RAN2)
C-XXX-C
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-C
C-XXX-C     Numerical Recipes Program For Generating Random Numbers.
C-XXX-C     The Common Block Parameters Are Used To Keep Track Of Info
C-XXX-C     as if it were a SAVE statement.
C-XXX-c     caution: -iff- is undefined before first call?
C-XXX-c     caution: should there be integer*4 typing here?
C-XXX-c
C-XXX-      PARAMETER (M1=259200,IA1=7141,IC1=54773)
C-XXX-      PARAMETER (M2=134456,IA2=8121,IC2=28411)
C-XXX-      PARAMETER (M3=243000,IA3=4561,IC3=51349)
C-XXX-      parameter (RM1=3.8580247D-06,RM2=7.4373773D-06)
C-XXX-c
C-XXX-      common/rands/ IFF,IX1,IX2,IX3,R(97)
C-XXX-c
C-XXX-      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
C-XXX-        IFF=1
C-XXX-        IX1=MOD(IC1-IDUM,M1)
C-XXX-        IX1=MOD(IA1*IX1+IC1,M1)
C-XXX-        IX2=MOD(IX1,M2)
C-XXX-        IX1=MOD(IA1*IX1+IC1,M1)
C-XXX-        IX3=MOD(IX1,M3)
C-XXX-        DO 11 J=1,97
C-XXX-          IX1=MOD(IA1*IX1+IC1,M1)
C-XXX-          IX2=MOD(IA2*IX2+IC2,M2)
C-XXX-          R(J)=(IX1+IX2*RM2)*RM1
C-XXX-   11   CONTINUE
C-XXX-        IDUM=1
C-XXX-      ENDIF
C-XXX-c
C-XXX-      IX1=MOD(IA1*IX1+IC1,M1)
C-XXX-      IX2=MOD(IA2*IX2+IC2,M2)
C-XXX-      IX3=MOD(IA3*IX3+IC3,M3)
C-XXX-      J=1+(97*IX3)/M3
C-XXX-      IF(J.GT.97.OR.J.LT.1) call abrt
C-XXX-      RAN2=R(J)
C-XXX-      R(J)=(IX1+IX2*RM2)*RM1
C-XXX-      RETURN
C-XXX-      END
C
C-XXX-c*module mdinit   *deck clrqtd
C-XXX-      subroutine clrqtd
C-XXX-c
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-c
C-XXX-      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-      parameter (zero=0.0d+00)
C-XXX-c
C-XXX-      common /atmdat/ vx(mxfrg), vy(mxfrg), vz(mxfrg),
C-XXX-     *                qw(mxfrg), qx(mxfrg), qy(mxfrg), qz(mxfrg),
C-XXX-     *                ox(mxfrg), oy(mxfrg), oz(mxfrg), qw1(mxfrg),
C-XXX-     *                qx1(mxfrg),qy1(mxfrg),qz1(mxfrg),qw2(mxfrg),
C-XXX-     *                qx2(mxfrg),qy2(mxfrg),qz2(mxfrg)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-C
C-XXX-c        **************************************************
C-XXX-c        *** clear quat derivative (1st and 2nd) arrays ***
C-XXX-c        **************************************************
C-XXX-c
C-XXX-       do i = 1, nfrg
C-XXX-         qw1(i) = zero
C-XXX-         qx1(i) = zero
C-XXX-         qy1(i) = zero
C-XXX-         qz1(i) = zero
C-XXX-         qw2(i) = zero
C-XXX-         qx2(i) = zero
C-XXX-         qy2(i) = zero
C-XXX-         qz2(i) = zero
C-XXX-       enddo
C-XXX-       return
C-XXX-       end
C
C-XXX-c*module mdinit   *deck pchk(ptol)
C-XXX-      subroutine pchk(ptol)
C-XXX-c        *********************************************************
C-XXX-c        *** momentum adjustment  ***
C-XXX-c        *********************************************************
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-C
C-XXX-      double precision MASTOT
C-XXX-      PARAMETER (MXPT=2000, MXFGPT=12000,MXFRG=1050,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-      parameter (zero=0.0d+00)
C-XXX-      dimension genmas(MXDFG)
C-XXX-c
C-XXX-      common /atmdat/ vx(mxfrg), vy(mxfrg), vz(mxfrg),
C-XXX-     *                qw(mxfrg), qx(mxfrg), qy(mxfrg), qz(mxfrg),
C-XXX-     *                ox(mxfrg), oy(mxfrg), oz(mxfrg), qw1(mxfrg),
C-XXX-     *                qx1(mxfrg),qy1(mxfrg),qz1(mxfrg),qw2(mxfrg),
C-XXX-     *                qx2(mxfrg),qy2(mxfrg),qz2(mxfrg)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
C-XXX-     *                FPNUC(MXPT),FGNUC(MXFGPT)
C-XXX-      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
C-XXX-c
C-XXX-      px=zero
C-XXX-      py=zero
C-XXX-      pz=zero
C-XXX-C
C-XXX-C mastot is total sum of all fmass(i)
C-XXX-      mastot=0
C-XXX-      do k=1,nfrg
C-XXX-         genmas(ISET(k)) = zero
C-XXX-         do i = 1,NMPTS(k)
C-XXX-            genmas(ISET(k)) = genmas(ISET(k)) + fmass(i)
C-XXX-         enddo
C-XXX-         mastot=mastot+genmas(ISET(k))
C-XXX-      enddo
C-XXX-C
C-XXX-      do 10 j=1,nfrg
C-XXX-         px=px+vx(j)*genmas(ISET(j))
C-XXX-         py=py+vy(j)*genmas(ISET(j))
C-XXX-         pz=pz+vz(j)*genmas(ISET(j))
C-XXX-   10 continue
C-XXX-      px=px/mastot
C-XXX-      py=py/mastot
C-XXX-      pz=pz/mastot
C-XXX-c        next silly line suppresses FTNCHEK warning on PTOL arg
C-XXX-      if(ptol.eq.1.234d+56) ptol=0.0d+00
C-XXX-c      if(abs(px).gt.ptol) goto 25
C-XXX-c      if(abs(py).gt.ptol) goto 25
C-XXX-c      if(abs(pz).gt.ptol) goto 25
C-XXX-c     return
C-XXX-c
C-XXX-      do 30 i=1,nfrg
C-XXX-        vx(i)=vx(i)-px
C-XXX-        vy(i)=vy(i)-py
C-XXX-        vz(i)=vz(i)-pz
C-XXX-   30 continue
C-XXX-      return
C-XXX-      end
C
C-XXX-c*module mdinit   *deck angvel
C-XXX-      subroutine angvel
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-C
C-XXX-      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-      parameter (two=2.0d+00)
C-XXX-c
C-XXX-      common /atmdat/ vx(mxfrg), vy(mxfrg), vz(mxfrg),
C-XXX-     *                qw(mxfrg), qx(mxfrg), qy(mxfrg), qz(mxfrg),
C-XXX-     *                ox(mxfrg), oy(mxfrg), oz(mxfrg), qw1(mxfrg),
C-XXX-     *                qx1(mxfrg),qy1(mxfrg),qz1(mxfrg),qw2(mxfrg),
C-XXX-     *                qx2(mxfrg),qy2(mxfrg),qz2(mxfrg)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      common /mdvel / wx(mxfrg),wy(mxfrg),wz(mxfrg),tqi(3,mxfrg)
C-XXX-c
C-XXX-c        ***********************************************
C-XXX-c        *** find initial ang. vel., w(0) from quat. ***
C-XXX-c        ***********************************************
C-XXX-c
C-XXX-      call vclr(wx,1,mxfrg)
C-XXX-      call vclr(wy,1,mxfrg)
C-XXX-      call vclr(wz,1,mxfrg)
C-XXX-c
C-XXX-      do 10 i=1,nfrg
C-XXX-       wx(i)=two*(-qz(i)*qy1(i)+qw(i)*qx1(i)+qy(i)*qz1(i)-
C-XXX-     *         qx(i)*qw1(i))
C-XXX-       wy(i)=two*(-qw(i)*qy1(i)-qz(i)*qx1(i)+qx(i)*qz1(i)-
C-XXX-     *         qy(i)*qw1(i))
C-XXX-       wz(i)=two*(qx(i)*qy1(i)-qy(i)*qx1(i)+qw(i)*qz1(i)-
C-XXX-     *         qz(i)*qw1(i))
C-XXX-   10 continue
C-XXX-      return
C-XXX-      end
C
C*MODULE MDINIT   *DECK MOVEVV(INTMD)
      SUBROUTINE MOVEVV(INTMD)
C       *********************************************************
C       *** Integrate equations of motion using a velocity    ***
C       *** verlet algorithm and quaternions                  ***
C       *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MDTYP,IRATYZ,IRATZX,IRATXY
C
      PARAMETER (MXFGPT=12000, MXFRG=1050,
     *           MXPT=2000, MXDFG=5, MXDPPT=MXFRG*MXDFG*12,
     *           MXAO=8192, MXSH=5000, MXGTOT=20000)
C
      PARAMETER (HALF=0.5D+00,TWO=2.0D+00)
C
      DIMENSION WX0(MXFRG),WY0(MXFRG),WZ0(MXFRG)
      DIMENSION IRATYZ(MXDFG),IRATZX(MXDFG),IRATXY(MXDFG)
      DIMENSION CL(MXFRG),PX(MXFRG),PY(MXFRG),PZ(MXFRG),QQQ(MXFRG)
      DIMENSION PWB(MXFRG),PXB(MXFRG),PYB(MXFRG),PZB(MXFRG)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FMCOM / XX(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C
      EKINT = 0.0D+00
      EKINR = 0.0D+00
      NMST = 0
      NPST = 0
      NDPST = 0
      NRST = 0
      NATST = 0
      NORST = 0
      NDST = 0
      NAOST = 0
C
      MXBF=0
C
      DO I=1,NFRG
         MXBF=MAX(MXBF,NPBF(I))
      END DO
      CALL VALFM(LOADFM)
      LPROVEC = LOADFM  + 1
      LCANVEC = LPROVEC + MXBF*NTMO
      LAST    = LCANVEC + MXBF*NTAO
C
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
      LENPV=MXBF*NTMO
      IF (NFRG.GT.0.AND.NTMO.GT.0)
     *   CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
      LENCV=MXBF*NTAO
      IF (NFRG.GT.0.AND.NTAO.GT.0)
     *   CALL DAREAD(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
C
      DT6=DT*DT/TWO
      DT2=DT/TWO
C
C
C      Translational and Rotational KE at t
C       do 110 i=1,nfrg
C        ekint = ekint + vx(i)**2 + vy(i)**2 + vz(i)**2
C        ekinr = ekinr + amom(1) * wx(i)**2 + amom(2) * wy(i)**2
C     *                + amom(3) * wz(i)**2
C  110 continue
C
C       ekint = half * ek2kcal * totmas * ekint
C       ekinr = half * ek2kcal * ekinr
C
           DO 4001 J=1,NFRG
C
C
      IRATYZ(ISET(J))=(AMOM(2,ISET(J))-
     *                AMOM(3,ISET(J)))/AMOM(1,ISET(J))
      IRATZX(ISET(J))=(AMOM(3,ISET(J))-
     *                 AMOM(1,ISET(J)))/AMOM(2,ISET(J))
      IRATXY(ISET(J))=(AMOM(1,ISET(J))-
     *                 AMOM(2,ISET(J)))/AMOM(3,ISET(J))
C       generate the correction coefficients
           S1=QW1(J)**2+QX1(J)**2+QY1(J)**2+QZ1(J)**2
           S2=QW1(J)*QW2(J)+QX1(J)*QX2(J)+QY1(J)*QY2(J)+QZ1(J)*QZ2(J)
           S3=QW2(J)**2+QX2(J)**2+QY2(J)**2+QZ2(J)**2
           S4=1-DT*DT*S1-DT*DT*DT*S2-DT6*DT6*(S3-S1*S1)
           CL(J)=1.0D+00-DT6*S1-SQRT(S4)
4001    CONTINUE
C
C
      DO 40 J=1,NFRG
C     Move positions to t+dt. Note
C     r(t+dt) = r(t)+dt*v(t)+1/2*dt^2*a(t)
C     AND r(t) is included in mddspq.src
        CENTDX=(VX(J)*DT+DT6*DEFT(1,J)*TOTMAI(ISET(J)))
        CENTDY=(VY(J)*DT+DT6*DEFT(2,J)*TOTMAI(ISET(J)))
        CENTDZ=(VZ(J)*DT+DT6*DEFT(3,J)*TOTMAI(ISET(J)))
C
        VX(J)=VX(J)+DT2*DEFT(1,J)*TOTMAI(ISET(J))
        VY(J)=VY(J)+DT2*DEFT(2,J)*TOTMAI(ISET(J))
        VZ(J)=VZ(J)+DT2*DEFT(3,J)*TOTMAI(ISET(J))
C
C   Save quaternions at time t
        QWSAV = QW(J)
        QXSAV = QX(J)
        QYSAV = QY(J)
        QZSAV = QZ(J)
C      write(6,*) 'Frag ',j,' Old Q: ',qwsav,qxsav,qysav,qzsav
C     Integrate quaternions, form Q(t+dt), add corrections:
C     Q(t+dt) = Q(t) + Q1d(t)*dt + 1/2*dt^2*Q2d(t) - c1*Q(t)
       QW(J)=QW(J)+DT*QW1(J)+DT6*QW2(J)-CL(J)*QW(J)
       QX(J)=QX(J)+DT*QX1(J)+DT6*QX2(J)-CL(J)*QX(J)
       QY(J)=QY(J)+DT*QY1(J)+DT6*QY2(J)-CL(J)*QY(J)
       QZ(J)=QZ(J)+DT*QZ1(J)+DT6*QZ2(J)-CL(J)*QZ(J)
C       write(6,*) 'Frag ',j,' New Q: ',qw(j),qx(j),qy(j),qz(j)
C
C  Do rotation and translation to get ecp to t+dt
C  Note: we need to call efcm before coming back
C        here to get newest COM coords
C
       CALL MDDSPQ(J,CENTDX,CENTDY,CENTDZ,QWSAV,QXSAV,QYSAV,QZSAV,
     *             QW(J),QX(J),QY(J),QZ(J),
     *             NMST,NPST,NDPST,NRST,NDST,NATST,NORST,XX(LPROVEC),
     *             NAOST,XX(LCANVEC),MXBF,INTMD)
C
      NMST  = NMST  + NMPTS(J)
      NPST  = NPST  + NPPTS(J)
      NDPST = NDPST + NDPPTS(J)
      NRST  = NRST  + NRPTS(J)
      NDST  = NDST  + NDPTS(J)
      NATST = NATST + NATEF(ISET(J))
      NORST = NORST + NORB(J)
      NAOST = NAOST + NAO(ISET(J))
C
   40 CONTINUE
C
C
      IF (NFRG.GT.0.AND.NTMO.GT.0)
     *     CALL DAWRIT(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
      IF (NFRG.GT.0.AND.NTAO.GT.0)
     *     CALL DAWRIT(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
      CALL RETFM(NEED)
C
C     Form inertia ratio, part of zeroth est. of w0:
C     w0 = w(t) + 1/2*dt*w1d(t)
C     where w1d(t) = torq/imom + irat*w*w
      DO 60 I=1,NFRG
       PX(I)=IRATYZ(ISET(I))*WY(I)*WZ(I)
       PY(I)=IRATZX(ISET(I))*WZ(I)*WX(I)
       PZ(I)=IRATXY(ISET(I))*WX(I)*WY(I)
       WX0(I)=WX(I)+DT2*(TQI(1,I)+PX(I))
       WY0(I)=WY(I)+DT2*(TQI(2,I)+PY(I))
       WZ0(I)=WZ(I)+DT2*(TQI(3,I)+PZ(I))
   60 CONTINUE
C
C        Calculate energy and forces at new t+dt
C
C       if (mdtyp.eq.efp) then
           CALL ENEDEF
           CALL TORQMO
C       else
C         call enedab
C       endif
C
       DO 80 I=1,NFRG
C
C      Move velocity from t+1/2dt to t+dt:
C      v(t+dt) = v(t+1/2dt) + 1/2*dt*a(t+dt)
        VX(I) = VX(I)+DT2*DEFT(1,I)*TOTMAI(ISET(I))
        VY(I) = VY(I)+DT2*DEFT(2,I)*TOTMAI(ISET(I))
        VZ(I) = VZ(I)+DT2*DEFT(3,I)*TOTMAI(ISET(I))
C
C      Change torq(t+dt) --> torq(t+dt)/amom
C        tqi(1,i)=torq(1,i)/amom(1)
C        tqi(2,i)=torq(2,i)/amom(2)
C        tqi(3,i)=torq(3,i)/amom(3)
C
C      Form first part of quat 2nd derivative:
C      Qd2(t+dt) = 1/2 * M * (torq/amom)
         QW2(I)=0.5D+00*(-QY(I)*TQI(1,I)-QZ(I)*TQI(2,I)+
     1         QX(I)*TQI(3,I))
         QX2(I)=0.5D+00*( QZ(I)*TQI(1,I)-QY(I)*TQI(2,I)-
     2         QW(I)*TQI(3,I))
         QY2(I)=0.5D+00*(QW(I)*TQI(1,I)+QX(I)*TQI(2,I)+
     3         QZ(I)*TQI(3,I))
         QZ2(I)=0.5D+00*(-QX(I)*TQI(1,I)+QW(I)*TQI(2,I)-
     4         QY(I)*TQI(3,I))
C
C      Form zeroth est. for w0 at t+dt:
C      w0(t+dt) = w0 + 1/2 * dt * tqi(t+dt)/amom
        WX0(I)=WX0(I)+DT2*TQI(1,I)
        WY0(I)=WY0(I)+DT2*TQI(2,I)
        WZ0(I)=WZ0(I)+DT2*TQI(3,I)
C
C      New px, py, pz with w0(t+dt) estimates
        PX(I)=IRATYZ(ISET(I))*WY0(I)*WZ0(I)
        PY(I)=IRATZX(ISET(I))*WZ0(I)*WX0(I)
        PZ(I)=IRATXY(ISET(I))*WX0(I)*WY0(I)
C
C      Form new angular velocity and p:
C      w(t+dt) = w0(t+dt) + 1/2 * dt * irat*w*w
        DO 350 L=1,3
           WX(I)=WX0(I)+PX(I)*DT2
           WY(I)=WY0(I)+PY(I)*DT2
           WZ(I)=WZ0(I)+PZ(I)*DT2
           PX(I)=IRATYZ(ISET(I))*WY(I)*WZ(I)
           PY(I)=IRATZX(ISET(I))*WZ(I)*WX(I)
           PZ(I)=IRATXY(ISET(I))*WX(I)*WY(I)
  350   CONTINUE
   80 CONTINUE
C
C
        DO 4005 I=1,NFRG
C      Generate the first time derivatives of the Q's.
         QW1(I)=(-QY(I)*WX(I)-QZ(I)*WY(I)+QX(I)*WZ(I) )*0.5D+00
         QX1(I)=( QZ(I)*WX(I)-QY(I)*WY(I)-QW(I)*WZ(I) )*0.5D+00
         QY1(I)=( QW(I)*WX(I)+QX(I)*WY(I)+QZ(I)*WZ(I) )*0.5D+00
         QZ1(I)=(-QX(I)*WX(I)+QW(I)*WY(I)-QY(I)*WZ(I) )*0.5D+00
          QQQ(I)=QW1(I)*QW1(I)+QX1(I)*QX1(I)+
     Q            QY1(I)*QY1(I)+QZ1(I)*QZ1(I)
         PX(I)=IRATYZ(ISET(I))*WY(I)*WZ(I)
         PY(I)=IRATZX(ISET(I))*WZ(I)*WX(I)
         PZ(I)=IRATXY(ISET(I))*WX(I)*WY(I)
            PWB(I)=0.5D+00*(-QY(I)*PX(I)-QZ(I)*PY(I)+
     1       QX(I)*PZ(I))
            PXB(I)=0.5D+00*( QZ(I)*PX(I)-QY(I)*PY(I)-
     2       QW(I)*PZ(I))
            PYB(I)=0.5D+00*( QW(I)*PX(I)+QX(I)*PY(I)+
     3       QZ(I)*PZ(I))
            PZB(I)=0.5D+00*(-QX(I)*PX(I)+QW(I)*PY(I)-
     4       QY(I)*PZ(I))
C        Complete the second time derivatives of the Q's
         QW2(I)=QW2(I)-QW(I)*QQQ(I)+PWB(I)
         QX2(I)=QX2(I)-QX(I)*QQQ(I)+PXB(I)
         QY2(I)=QY2(I)-QY(I)*QQQ(I)+PYB(I)
         QZ2(I)=QZ2(I)-QZ(I)*QQQ(I)+PZB(I)
C       qchk(i)=qw(i)**2+qx(i)**2+qy(i)**2+qz(i)**2
C       qchk2(i)=qw(i)*qw1(i)+qx(i)*qx1(i)+qy(i)*qy1(i)+
C     *            qz(i)*qz1(i)
C        if(i.le.5) then
C       write(6,*) 'i,qchk,qchk2 = ',i,qchk(i),qchk2(i)
C        endif
4005   CONTINUE
C
C
C      Translational and Rotational KE at t+dt
       DO 110 I=1,NFRG
        EKINT = EKINT + TOTMAS(ISET(I))*(VX(I)**2 + VY(I)**2 + VZ(I)**2)
        EKINR = EKINR + AMOM(1,ISET(I)) * WX(I)**2 + AMOM(2,ISET(I))
     *              * WY(I)**2 + AMOM(3,ISET(I)) * WZ(I)**2
  110 CONTINUE
C
C                ekint = half * ek2kcal * TOTMAS * ekint
       EKINT = HALF * EK2KCAL * EKINT
       EKINR = HALF * EK2KCAL * EKINR
       RETURN
       END
C*MODULE MDINIT   *DECK TORQMO
       SUBROUTINE TORQMO
C        *********************************************************
C        *** torque to torque/moment of inertia  ***
C        *********************************************************
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFGPT=12000,MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
C
C      Change torq(t+dt) --> torq(t+dt)/amom
C
       DO I=1,NFRG
         AXX=-QW(I)**2+QX(I)**2-QY(I)**2+QZ(I)**2
         AXY= 2.0D+00*(QY(I)*QZ(I)-QW(I)*QX(I))
         AXZ= 2.0D+00*(QX(I)*QY(I)+QW(I)*QZ(I))
C
         AYX=-2.0D+00*(QW(I)*QX(I)+QY(I)*QZ(I))
         AYY= QW(I)**2-QX(I)**2-QY(I)**2+QZ(I)**2
         AYZ= 2.0D+00*(QX(I)*QZ(I)-QW(I)*QY(I))
C
         AZX= 2.0D+00*(QX(I)*QY(I)-QW(I)*QZ(I))
         AZY=-2.0D+00*(QW(I)*QY(I)+QX(I)*QZ(I))
         AZZ=-QW(I)**2-QX(I)**2+QY(I)**2+QZ(I)**2
C
        TQI(1,I)=AXX*TORQ(1,I)+AXY*TORQ(2,I)+AXZ*TORQ(3,I)
        TQI(2,I)=AYX*TORQ(1,I)+AYY*TORQ(2,I)+AYZ*TORQ(3,I)
        TQI(3,I)=AZX*TORQ(1,I)+AZY*TORQ(2,I)+AZZ*TORQ(3,I)
C
C
        TQI(1,I)=TQI(1,I)/AMOM(1,ISET(I))
        TQI(2,I)=TQI(2,I)/AMOM(2,ISET(I))
        TQI(3,I)=TQI(3,I)/AMOM(3,ISET(I))
       ENDDO
      RETURN
      END
C
C-XXX-c*module mdinit   *deck gaussw
C-XXX-      SUBROUTINE GAUSSW(temp)
C-XXX-c
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-C
C-XXX-      PARAMETER (MXFGPT=12000,MXFRG=1050,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-C
C-XXX-      common /convmd/ au2kcal, ek2kcal, au2ang, au2sec, au2amu
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
C-XXX-      common /mdvel / wx(mxfrg),wy(mxfrg),wz(mxfrg),tqi(3,mxfrg)
C-XXX-      common /moddat/ amom(3,MXDFG), totmas(MXDFG), totmai(MXDFG),
C-XXX-     *                  drx(MXFGPT), dry(MXFGPT), drz(MXFGPT)
C-XXX-c
C-XXX-      DATA HALF,ONE,TWO/0.5d+00,1.0d+00,2.0d+00/
C-XXX-c
C-XXX-c         initialize RNG
C-XXX-c
C-XXX-      iseed = -43891265
C-XXX-      call ranth(iseed,rannum)
C-XXX-      iseed = 0
C-XXX-c
C-XXX-c      m2bohr=1.0d+00/5.2917706d-11
C-XXX-      boltzk = 0.83112120d+00/1.0d-12/1.0d-12/au2ang/au2ang
C-XXX-      TWOPI=TWO*ACOS(-ONE)
C-XXX-c
C-XXX-      DO 2 I=1,nfrg
C-XXX-          SIGMAX=SQRT((BOLTZK*TEMP)/amom(1,ISET(i)))
C-XXX-          SIGMAY=SQRT((BOLTZK*TEMP)/amom(2,ISET(i)))
C-XXX-          SIGMAZ=SQRT((BOLTZK*TEMP)/amom(3,ISET(i)))
C-XXX-          CALL RANth(ISEED,R1)
C-XXX-          CALL RANth(ISEED,R2)
C-XXX-          FACT=SQRT(-TWO*LOG(R1))
C-XXX-          RANN1=FACT*COS(TWOPI*R2)
C-XXX-          RANN2=FACT*SIN(TWOPI*R2)
C-XXX-          CALL RANth(ISEED,R1)
C-XXX-          CALL RANth(ISEED,R2)
C-XXX-          FACT=SQRT(-TWO*LOG(R1))
C-XXX-          RANN3=FACT*COS(TWOPI*R2)
C-XXX-          CALL RANth(ISEED,R1)
C-XXX-          CALL RANth(ISEED,R2)
C-XXX-          CALL RANth(ISEED,R3)
C-XXX-          IF (R1.GT.HALF)RANN1=-RANN1
C-XXX-          IF (R2.GT.HALF)RANN2=-RANN2
C-XXX-          IF (R3.GT.HALF)RANN3=-RANN3
C-XXX-          WX(I)=SIGMAX*RANN1
C-XXX-          WY(I)=SIGMAY*RANN2
C-XXX-          WZ(I)=SIGMAZ*RANN3
C-XXX-C    1   CONTINUE
C-XXX-    2 CONTINUE
C-XXX-      RETURN
C-XXX-      END
C
C*MODULE MDINIT   *DECK RANBS(ISEED)
      SUBROUTINE RANBS(ISEED,RANDOM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER*4 I,I1,I2,I3,I4,I5,I6,I7,I8,I9
      DATA I2/16807/,I4/32768/,I5/65536/,I3/2147483647/
C
C           algorithm is based on 4-byte values, but the
C           seed's data type is left to compiler flags.
C
      I=ISEED
C
      I6=I/I5
      I7=(I-I6*I5)*I2
      I8=I7/I5
      I9=I6*I2+I8
      I1=I9/I4
      I=(((I7-I8*I5)-I3)+(I9-I1*I4)*I5)+I1
      IF(I.LT.0) I=I+I3
C
      RANDOM=(4.656613D-10)*I
      ISEED=I
      RETURN
      END
C
C*MODULE MDINIT   *DECK NORMT
      SUBROUTINE NORMT(TEMP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFGPT=12000,MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION SET(MXFRG*3)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
      SAVE ISEED
C
      BOLTZK = 0.83112120D+00/1.0D-12/1.0D-12/AU2ANG/AU2ANG
C
C     for u1,u2 uniformly distrib on (0,1)
C     x = sqrt(-2*ln(u1)) * cos(2*pi*u2)
C     gives normally distrib x on (-inf,+inf)
C
      ISEED=1234567
C
      PI=4.0D+00*ATAN(1.0D+00)
      NTEST=NFRG*3
      XMEAN=0.0D+00
C
      DO I=1,NTEST
        CALL RANBS(ISEED,U1)
        CALL RANBS(ISEED,U2)
        SET(I) = SQRT(-2.0D+00*LOG(U1)) * COS(2.0D+00*PI*U2)
        XMEAN=XMEAN+SET(I)
      ENDDO
      XMEAN=XMEAN/NTEST
C
      XMSD =0.0D+00
      DO I=1,NTEST
        XMSD=XMSD+(SET(I)-XMEAN)**2
      ENDDO
      XMSD=XMSD/NTEST
C
C      write(*,'('' Mean and RMSD = '',2f10.5)') xmean,sqrt(xmsd)
C
C     we multiply the random number by sqrt(kT/m) to
C     get a velocity in units of sqrt(kcal/(mole-amu))
C     this is our internal velocity unit so that
C     sum over m*v**2 has units of kcal/mole
C
C
C      energy=0.0d+00
      J=0
      DO I=1,NFRG
        SIGMA=SQRT(BOLTZK*TEMP/FRGMAS(I))
        VX(I)=SIGMA*SET(J+1)
        VY(I)=SIGMA*SET(J+2)
        VZ(I)=SIGMA*SET(J+3)
        J=J+3
C        energy=energy+0.5d+00*hmass*veloc**2
      ENDDO
C
C     compute temperature from this KE
C     corresponds to ntest degrees of freedom
C      temp = energy/(0.5d+00*boltz*ntest)
C      write(*,'('' Temperature = '',f10.5)') temp
C
C
      RETURN
      END
C
C*MODULE MDINIT   *DECK NORMW
      SUBROUTINE NORMW(TEMP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFGPT=12000,MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION SET(MXFRG*3)
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                  DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
C
      SAVE ISEED
C
      BOLTZK = 0.83112120D+00/1.0D-12/1.0D-12/AU2ANG/AU2ANG
C
C     for u1,u2 uniformly distrib on (0,1)
C     x = sqrt(-2*ln(u1)) * cos(2*pi*u2)
C     gives normally distrib x on (-inf,+inf)
C
      ISEED=7654321
C
      PI=4.0D+00*ATAN(1.0D+00)
      NTEST=NFRG*3
      XMEAN=0.0D+00
C
      DO I=1,NTEST
        CALL RANBS(ISEED,U1)
        CALL RANBS(ISEED,U2)
        SET(I) = SQRT(-2.0D+00*LOG(U1)) * COS(2.0D+00*PI*U2)
        XMEAN=XMEAN+SET(I)
      ENDDO
      XMEAN=XMEAN/NTEST
C
      XMSD =0.0D+00
      DO I=1,NTEST
        XMSD=XMSD+(SET(I)-XMEAN)**2
      ENDDO
      XMSD=XMSD/NTEST
C
C      write(*,'('' Mean and RMSD = '',2f10.5)') xmean,sqrt(xmsd)
C
C     we multiply the random number by sqrt(kT/m) to
C     get a velocity in units of sqrt(kcal/(mole-amu))
C     this is our internal velocity unit so that
C     sum over m*v**2 has units of kcal/mole
C
C
C      energy=0.0d+00
      J=0
      DO I=1,NFRG
         SIGMAX=SQRT(BOLTZK*TEMP/AMOM(1,ISET(I)))
         SIGMAY=SQRT(BOLTZK*TEMP/AMOM(2,ISET(I)))
         SIGMAZ=SQRT(BOLTZK*TEMP/AMOM(3,ISET(I)))
         WX(I)=SIGMAX*SET(J+1)
         WY(I)=SIGMAY*SET(J+2)
         WZ(I)=SIGMAZ*SET(J+3)
         J=J+3
C        energy=energy+0.5d+00*hmass*veloc**2
      ENDDO
C
C     compute temperature from this KE
C     corresponds to ntest degrees of freedom
C      temp = energy/(0.5d+00*boltz*ntest)
C      write(*,'('' Temperature = '',f10.5)') temp
C
C
      RETURN
      END
C
C*MODULE MDINIT   *DECK EU2QUT
      SUBROUTINE EU2QUT
C       *********************************************************
C       *** get correct quaternions from coordinates and    *****
C       *** Euler angles.                                     ***
C       *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFGPT=12000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      CHARACTER*8 FRGNME
      DIMENSION R(3,3),VSCRATCH(9),COORD(3,3)
C---  dimension rt(3,3)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
C
      IC1 = 0
      DO I = 1,NFRG
         IC1 = IC1 + 1
C       first get rotation matrix from molecule frame to lab frame
C
C       Compute the Euler angles between lab frame and molecule frame
C
C       get unit vector parallel to dipole moment
         COORD(1,1) = EFC(1,IC1) - EFCENT(1,I)
        COORD(1,2) = EFC(2,IC1) - EFCENT(2,I)
        COORD(1,3) = EFC(3,IC1) - EFCENT(3,I)
        COORD(2,1) = EFC(1,IC1+1) - EFCENT(1,I)
        COORD(2,2) = EFC(2,IC1+1) - EFCENT(2,I)
        COORD(2,3) = EFC(3,IC1+1) - EFCENT(3,I)
        COORD(3,1) = EFC(1,IC1+2) - EFCENT(1,I)
        COORD(3,2) = EFC(2,IC1+2) - EFCENT(2,I)
        COORD(3,3) = EFC(3,IC1+2) - EFCENT(3,I)
C
        UX=COORD(2,1)+COORD(3,1)-2.0D+00*COORD(1,1)
        UY=COORD(2,2)+COORD(3,2)-2.0D+00*COORD(1,2)
        UZ=COORD(2,3)+COORD(3,3)-2.0D+00*COORD(1,3)
C
C
C
C        ux=efc(1,ic1+1)+efc(1,ic1+2)-2.0d+00*efc(1,ic1)
C        uy=efc(2,ic1+1)+efc(2,ic1+2)-2.0d+00*efc(2,ic1)
C        uz=efc(3,ic1+1)+efc(3,ic1+2)-2.0d+00*efc(3,ic1)
        U=SQRT(UX**2+UY**2+UZ**2)
        UX=-UX/U
        UY=-UY/U
        UZ=-UZ/U
C       vector (ux,uy,uz) is now in direction of principle +z axis
C       (this is opposite in direction to dipole moment)
C
C       theta is angle in range 0 to pi between principle axis
C       and +z lab axis,
C       cos(theta)=uz; sin(theta)=sqrt(ux**2+uy**2)
C
        THETA=ATAN2(SQRT(UX**2+UY**2),UZ)
C
C       phi is angle in range -pi to pi about the lab frame z axis
C       between the lab frame x axis and the projection of the molecule
C       frame z axis onto the xy lab frame plane
C---    phix=atan2(uy,ux)
C       phi is angle in range -pi to pi about the lab frame z axis
C       between the lab frame y axis and the projection of the molecule
C       frame z axis onto the xy lab frame plane
        PHIY=-ATAN2(UX,UY)
C
C       following computes chi by applying a pair of rotations designed
C       to orient the molecule dipole along the lab frame z axis
C       first rotate by -phiy about z; then rotate by theta about x
        PHI=PHIY
        R(1,1)= COS(PHI)
        R(2,1)=-COS(THETA)*SIN(PHI)
        R(3,1)=-SIN(THETA)*SIN(PHI)
        R(1,2)= SIN(PHI)
        R(2,2)= COS(THETA)*COS(PHI)
        R(3,2)= SIN(THETA)*COS(PHI)
        R(1,3)= 0.0D+00
        R(2,3)=-SIN(THETA)
        R(3,3)= COS(THETA)
C       translate the molecule so its COM is at origin,
C       actually, the data in coord is ALREADY reletive to COM
C       apply rotation matrix
        VSCRATCH(1)=COORD(1,1)
        VSCRATCH(2)=COORD(1,2)
        VSCRATCH(3)=COORD(1,3)
        VSCRATCH(4)=COORD(2,1)
        VSCRATCH(5)=COORD(2,2)
        VSCRATCH(6)=COORD(2,3)
        VSCRATCH(7)=COORD(3,1)
        VSCRATCH(8)=COORD(3,2)
        VSCRATCH(9)=COORD(3,3)
C
C
        CALL ROTBILL(R,VSCRATCH)
C
C       now determine the chi angle to rotate about the z axis
C       to put H2 in the y-z plane of the molecule frame
C
        CHI=-ATAN2( VSCRATCH(7),VSCRATCH(8))
C       following matrix rotates by -chi about +z
        R(1,1)= COS(CHI)
        R(2,1)=-SIN(CHI)
        R(3,1)= 0.0D+00
        R(1,2)= SIN(CHI)
        R(2,2)= COS(CHI)
        R(3,2)= 0.0D+00
        R(1,3)= 0.0D+00
        R(2,3)= 0.0D+00
        R(3,3)= 1.0D+00
        CALL ROTBILL(R,VSCRATCH)
C
C       using the computed theta, phi, chi angles, construct
C       the matrix that rotates from the molecule frame to lab frame
C       the matrix that changes a vector from the lab basis to the
C       principle axis basis
C
C        Use zxz convention
          PHI=PHIY
          THETA=-THETA
C         following matrix rotates by phi about +z axis
C         then by theta about the new +x axis
C         then by chi about the new +z axis
          R(1,1)=-COS(THETA)*SIN(PHI)*SIN(CHI)+COS(PHI)*COS(CHI)
          R(2,1)=-COS(THETA)*SIN(PHI)*COS(CHI)-COS(PHI)*SIN(CHI)
          R(3,1)= SIN(THETA)*SIN(PHI)
          R(1,2)=+COS(THETA)*COS(PHI)*SIN(CHI)+SIN(PHI)*COS(CHI)
          R(2,2)=+COS(THETA)*COS(PHI)*COS(CHI)-SIN(PHI)*SIN(CHI)
          R(3,2)=-SIN(THETA)*COS(PHI)
          R(1,3)= SIN(THETA)*SIN(CHI)
          R(2,3)= SIN(THETA)*COS(CHI)
          R(3,3)= COS(THETA)
C
C       the rotation matrices constructed should transform the lab
C       frame displacement vectors into the principle molecule frame
C
C       use this matrix to rotate from molecule frame to lab frame
C---    do j=1,3
C---      do k=1,3
C---         rt(j,k)=r(k,j)
C---      enddo
C---    enddo
C
C       evaluate quaternions
C
C       following relationships are in Mountain's and Evans' papers
C       CONFIRMED that the Euler angles to use in these expressions
C       are those from the zxz convention.
        QCHI =COS(THETA/2.0D+00)*COS((CHI+PHI)/2.0D+00)
        QETA =SIN(THETA/2.0D+00)*COS((CHI-PHI)/2.0D+00)
        QXI  =SIN(THETA/2.0D+00)*SIN((CHI-PHI)/2.0D+00)
        QZETA=COS(THETA/2.0D+00)*SIN((CHI+PHI)/2.0D+00)
C
C       following are in Mountain's and Evan's papers
        QW(I)=QXI
        QX(I)=QETA
        QY(I)=QZETA
        QZ(I)=QCHI
C
C        if(i.le.6) then
C       following is from Evans
C        rq(1,1) =-qxi**2+qeta**2-qzeta**2+qchi**2
C        rq(2,1) = -2*(qxi*qeta+qzeta*qchi)
C        rq(3,1) =  2*(qeta*qzeta-qxi*qchi)
C
C        rq(1,2) =  2*(qzeta*qchi-qxi*qeta)
C        rq(2,2) =+qxi**2-qeta**2-qzeta**2+qchi**2
C        rq(3,2) = -2*(qxi*qzeta+qeta*qchi)
C
C        rq(1,3) =  2*(qeta*qzeta+qxi*qchi)
C        rq(2,3) =  2*(qeta*qchi-qxi*qzeta)
C        rq(3,3) =-qxi**2-qeta**2+qzeta**2+qchi**2
C
        IC1 = IC1+NMPTS(I)-1
        ENDDO
C
        RETURN
        END
C
C*MODULE MDINIT   *DECK ROTBILL(R,C)
      SUBROUTINE ROTBILL(R,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION R(3,3),C(9),CDUM(3)
C
C        ***********************************************************
C        *** utility subroutine to apply 3-by-3 rotation matrix  ***
C        ***  to three vectors                                   ***
C        ***********************************************************
C
      DO I=1,3
        SUM=0.0D+00
        DO J=1,3
          SUM=SUM+R(I,J)*C(J)
        ENDDO
        CDUM(I)=SUM
      ENDDO
      DO I=1,3
        C(I)=CDUM(I)
      ENDDO
C
      DO I=1,3
        SUM=0.0D+00
        DO J=1,3
          SUM=SUM+R(I,J)*C(3+J)
        ENDDO
        CDUM(I)=SUM
      ENDDO
      DO I=1,3
        C(3+I)=CDUM(I)
      ENDDO
C
      DO I=1,3
        SUM=0.0D+00
        DO J=1,3
          SUM=SUM+R(I,J)*C(6+J)
        ENDDO
        CDUM(I)=SUM
      ENDDO
      DO I=1,3
        C(6+I)=CDUM(I)
      ENDDO
C
      RETURN
      END
C
C*MODULE MDINIT   *DECK INQDOT
      SUBROUTINE INQDOT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION IRATYZ,IRATZX,IRATXY
C
      PARAMETER (MXFGPT=12000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION PX(MXFRG),PY(MXFRG),PZ(MXFRG),QQQ(MXFRG)
      DIMENSION PWB(MXFRG),PXB(MXFRG),PYB(MXFRG),PZB(MXFRG)
      DIMENSION IRATYZ(MXDFG),IRATZX(MXDFG),IRATXY(MXDFG)
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /MODDAT/ AMOM(3,MXDFG), TOTMAS(MXDFG), TOTMAI(MXDFG),
     *                DRX(MXFGPT), DRY(MXFGPT), DRZ(MXFGPT)
C
C       *********************************************************
C       *** get initial qdot and qdot2 from quats calculated  ***
C       *** from euler angles.                                ***
C       *********************************************************
C
      CALL TORQMO
C
      DO I = 1,NFRG
C
      IRATYZ(ISET(I))=(AMOM(2,ISET(I))-
     *                 AMOM(3,ISET(I)))/AMOM(1,ISET(I))
      IRATZX(ISET(I))=(AMOM(3,ISET(I))-
     *                 AMOM(1,ISET(I)))/AMOM(2,ISET(I))
      IRATXY(ISET(I))=(AMOM(1,ISET(I))-
     *                 AMOM(2,ISET(I)))/AMOM(3,ISET(I))
C
C      Generate the first time derivatives of the Q's.
      QW1(I)=(-QY(I)*WX(I)-QZ(I)*WY(I)+QX(I)*WZ(I) )*0.5D+00
         QX1(I)=( QZ(I)*WX(I)-QY(I)*WY(I)-QW(I)*WZ(I) )*0.5D+00
         QY1(I)=( QW(I)*WX(I)+QX(I)*WY(I)+QZ(I)*WZ(I) )*0.5D+00
         QZ1(I)=(-QX(I)*WX(I)+QW(I)*WY(I)-QY(I)*WZ(I) )*0.5D+00
C
         QQQ(I)=QW1(I)*QW1(I)+QX1(I)*QX1(I)+
     *            QY1(I)*QY1(I)+QZ1(I)*QZ1(I)
C
         PX(I)=IRATYZ(ISET(I))*WY(I)*WZ(I)
         PY(I)=IRATZX(ISET(I))*WZ(I)*WX(I)
         PZ(I)=IRATXY(ISET(I))*WX(I)*WY(I)
         PWB(I)=0.5D+00*(-QY(I)*PX(I)-QZ(I)*PY(I)+
     1       QX(I)*PZ(I))
         PXB(I)=0.5D+00*( QZ(I)*PX(I)-QY(I)*PY(I)-
     2       QW(I)*PZ(I))
         PYB(I)=0.5D+00*( QW(I)*PX(I)+QX(I)*PY(I)+
     3       QZ(I)*PZ(I))
         PZB(I)=0.5D+00*(-QX(I)*PX(I)+QW(I)*PY(I)-
     4       QY(I)*PZ(I))
C
C      Form first part of quat 2nd derivative:
C      Qd2(t+dt) = 1/2 * M * (torq/amom)
        QW2(I)=0.5D+00*(-QY(I)*TQI(1,I)-QZ(I)*TQI(2,I)+
     1         QX(I)*TQI(3,I))
         QX2(I)=0.5D+00*( QZ(I)*TQI(1,I)-QY(I)*TQI(2,I)-
     2         QW(I)*TQI(3,I))
          QY2(I)=0.5D+00*(QW(I)*TQI(1,I)+QX(I)*TQI(2,I)+
     3         QZ(I)*TQI(3,I))
         QZ2(I)=0.5D+00*(-QX(I)*TQI(1,I)+QW(I)*TQI(2,I)-
     4         QY(I)*TQI(3,I))
C
C        Complete the second time derivatives of the Q's
         QW2(I)=QW2(I)-QW(I)*QQQ(I)+PWB(I)
         QX2(I)=QX2(I)-QX(I)*QQQ(I)+PXB(I)
         QY2(I)=QY2(I)-QY(I)*QQQ(I)+PYB(I)
         QZ2(I)=QZ2(I)-QZ(I)*QQQ(I)+PZB(I)
C
C        write(6,*)'in inqdot: i,pwb(i),tqi(1,i)=',i,pwb(i),tqi(1,i)
        ENDDO
C
        RETURN
        END
C
C*MODULE MDDSPQ   *DECK MDDSPQ
      SUBROUTINE MDDSPQ(N,XDEL,YDEL,ZDEL,QWOLD,QXOLD,QYOLD,QZOLD,
     *                  QW,QX,QY,QZ,NMST,NPST,NDPST,NRST,NDST,NATST,
     *                  NORST,PROVEC,NAOST,CANVEC,MXBF,INTMD)
C
C    ****************************************************************
C    ** modification of efdspl to take quaternions for generating  **
C    ** rotation matrix         (JMS)                              **
C    ****************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      PARAMETER (MXPT=2000, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXAO=8192, MXSH=5000, MXGTOT=20000,
     *           MXDPPT=MXFRG*MXDFG*12, MXIFRQ=12)
C
      PARAMETER (SQRT3=1.732050808D+00,SQRT5=2.236067978D+00)
      PARAMETER (SQRT7=2.645751311D+00)
      PARAMETER (TWO=2.0D+00,THREE=3.0D+00,FOUR=4.0D+00,SIX=6.0D+00)
      PARAMETER (TWELVE=12.0D+00)
      PARAMETER (CONST=2.0D+00/1.7320508D+00)
C
      CHARACTER*8 FRGNME,POLNAM,REPNAM,DPOLNAM
      CHARACTER*8 DNAME,DPSNAM
C
      DIMENSION T3(3,3),T3T(3,3)
      DIMENSION T1(3,3), T2(3,3)
      DIMENSION TEMP(3,MXPT),TEMP2(6),TEMP3(3),TEMP4(6,3),TEMP5(6,3)
      DIMENSION DTEMP(3,MXIFRQ*MXPT)
      DIMENSION TPXX(6),TPXY(6),TPYY(6),TPXZ(6),TPYZ(6),TPZZ(6)
      DIMENSION TPXX2(6),TPXY2(6),TPYY2(6),TPXZ2(6),TPYZ2(6),TPZZ2(6)
C
C      dimension prcopb(3,mxfgpt),centpb(3,mxfgpt)
      DOUBLE PRECISION MDTYP
      INTEGER INTMD
      LOGICAL PBCOUT,PBCWRI
C
      COMMON /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /EFPPAR/ EFP(3,MXFGPT),EFPOL(9,MXFGPT),
     *                EFDP(3,MXIFRQ*MXFGPT),EFDPOL(9,MXIFRQ*MXFGPT),
     *                ENO,DIND(3,MXFGPT),DINDD(3,MXFGPT),
     *                POLNAM(MXFGPT),DPOLNAM(MXIFRQ*MXFGPT)
      COMMON /EFPOTD/ METHOD,INABIO,MOVE
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /PBCCO/  FRGCPB(3,MXFGPT),EFPPB(3,MXFGPT),
     *                EFDPPB(3,MXFGPT),CREPPB(3,MXFGPT)
      COMMON /PBC/    PBCOUT,PBCWRI
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
      COMMON /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /REPPAR/ CREP(3,MXFGPT),CLPR(4*MXFGPT),ZLPR(4*MXFGPT),
     *                NLPR(4*MXFGPT),KFR(MXFGPT),KLR(MXFGPT),
     *                REPNAM(MXFGPT)
C
      DIMENSION PROVEC(MXBF,NTMO),CANVEC(MXBF,NTAO)
C
      NPTS = NMPTS(N)
C
C        write(6,*) 'in mddspq'
C
      CALL MDROTT(QWOLD,QXOLD,QYOLD,QZOLD,QW,QX,QY,QZ,T3,INTMD)
C
C      if(N.le.6) then
C      write(6,*) 'fragment ',N,'  UPDATE RTN MATRIX'
C      do i=1,3
C       write(6,*) t3(i,1), t3(i,2), t3(i,3)
C      enddo
C      endif
C
C      write(6,*) 't3(1,2,3) = ',t3(1,1),t3(2,1),t3(3,1)
C      CALL EFROTT(THETAX,THETAY,THETAZ,T3)
C
C  COMPUTE TRANSPOSE
C
      CALL TRPOSE(T3,T3T,3,3,0)
C
C   5.  ROTATE AND TRANSLATE STANDARD FRAGMENT TO COINCIDE WITH FRAG N
C
C  FIND FRAGMENT COORDINATES RELATIVE TO POINT OF ROTATION
C
C  find fragment coordinates relative to point of rotation
C  get BODY-FIXED coordinates--> COM = origin
C  Note that efcent and efc are at time t here...
C
      DO 505 I=1,NPTS
        DO 504 J=1,3
          TEMP(J,I) = EFC(J,I+NMST) - EFCENT(J,N)
  504   CONTINUE
  505 CONTINUE
C
C  NOW ROTATE FRAGMENT; rotate based on COM/body-fixed coordinates
C
      CALL MRARBR(T3,3,3,3,TEMP,3,NPTS,EFC(1,NMST+1),3)
C
C CHANGE COORDINATES RELATIVE TO EFCENT BACK TO OLD COORDINATE SYSTEM
C change back to SPACE-FIXED coordinates; note that efcent is at time t
C here
      DO 605 I=1,NPTS
        DO 604 J=1,3
          EFC(J,I+NMST) = EFC(J,I+NMST) + EFCENT(J,N)
  604   CONTINUE
  605 CONTINUE
C
C      write(6,*) 'xdel,y,z = ',xdel,ydel,zdel
C      write(6,*) 'efcent = ',efcent(1,n),efcent(2,n),efcent(3,n)
C  PBC??
      IF(PBCWRI) THEN
        XDELPB = XDEL - BOXL*ANINT(( XDEL + EFCENT(1,N) )*BOXLI)
        YDELPB = YDEL - BOXL*ANINT(( YDEL + EFCENT(2,N) )*BOXLI)
        ZDELPB = ZDEL - BOXL*ANINT(( ZDEL + EFCENT(3,N) )*BOXLI)
        IF(XDELPB.NE.XDEL.OR.YDELPB.NE.YDEL.OR.ZDELPB.NE.ZDEL) THEN
C          write(6,*) 'PBC!'
C          write(6,*) 'xdelpb,y,z =',xdelpb,ydelpb,zdelpb
        ENDIF
      ENDIF
C
C  FEED IN INFORMATION INTO FINAL ARRAYS AFTER ROTATION
C
C DO TRANSLATION
C so xdel is dt*v(t+1/2dt) from mdmove--we want to be moving the
C COM (since mdmove is in terms of COM motion)
C new efcent(t+dt) = efcent(t) + dt*v(t+1/2dt)
C new efc(t+dt) = efcent(t+dt) + rotated body-fixed atom positions,
C        AT*efc-bf(t)
C new efc(t+dt) = efcent(t) + dt*v(t+1/2dt) + rotated stuff
C We have added efcent(t) and rotated stuff above,
C so we just need to add in info from mdmove...
C
      IF (PBCWRI) THEN
      DO 449 I = 1,NPTS
        FRGCPB(1,NMST+I) = EFC(1,NMST+I) + XDELPB
        FRGCPB(2,NMST+I) = EFC(2,NMST+I) + YDELPB
        FRGCPB(3,NMST+I) = EFC(3,NMST+I) + ZDELPB
  449 CONTINUE
      ENDIF
C
C     WRITE(6,*)'ROTATED COORDINATES,ME,NSERCH',ME,NSERCH
      DO 550 I = 1,NPTS
        EFC(1,NMST+I) = EFC(1,NMST+I) + XDEL
        EFC(2,NMST+I) = EFC(2,NMST+I) + YDEL
        EFC(3,NMST+I) = EFC(3,NMST+I) + ZDEL
C       WRITE(6,FMT='(3F9.5)') (EFC(J,I+NMST), J=1,3)
  550 CONTINUE
      IF (METHOD .EQ. 6) RETURN
C
C ROTATING DIPOLE MOMENTS
C
      CALL MRARBR(T3,3,3,3,EFDIP(1,NMST+1),3,NPTS,TEMP,3)
      DO 555 I=1,NPTS
      DO 555 J=1,3
      EFDIP(J,NMST+I)=TEMP(J,I)
 555  CONTINUE
C
C ROTATING QUADRUPOLE MOMENTS
C
      TEMP3(1) = 0.0D+00
      TEMP3(2) = 0.0D+00
      TEMP3(3) = 0.0D+00
      DO 600 I=1,NPTS
        TEMP2(1)=EFQAD(1,I+NMST)
        TEMP2(2)=EFQAD(4,I+NMST)
        TEMP2(3)=EFQAD(2,I+NMST)
        TEMP2(4)=EFQAD(5,I+NMST)
        TEMP2(5)=EFQAD(6,I+NMST)
        TEMP2(6)=EFQAD(3,I+NMST)
        CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
        EFQAD(1,NMST+I)=TEMP4(1,1)
        EFQAD(4,NMST+I)=TEMP4(2,1)
        EFQAD(2,NMST+I)=TEMP4(3,1)
        EFQAD(5,NMST+I)=TEMP4(4,1)
        EFQAD(6,NMST+I)=TEMP4(5,1)
        EFQAD(3,NMST+I)=TEMP4(6,1)
  600 CONTINUE
C
C     Jay {{{
      IF(IEFP2.EQ.1) THEN
        CALL MRARBR(T3,3,3,3,DIPHF(1,NMST+1),3,NPTS,TEMP,3)
        DO 601 I=1,NPTS
          DO 601 J=1,3
            DIPHF(J,NMST+I)=TEMP(J,I)
 601    CONTINUE
C
        TEMP3(1) = 0.0D+00
        TEMP3(2) = 0.0D+00
        TEMP3(3) = 0.0D+00
        DO 602 I=1,NPTS
          TEMP2(1)=QADHF(1,I+NMST)
          TEMP2(2)=QADHF(4,I+NMST)
          TEMP2(3)=QADHF(2,I+NMST)
          TEMP2(4)=QADHF(5,I+NMST)
          TEMP2(5)=QADHF(6,I+NMST)
          TEMP2(6)=QADHF(3,I+NMST)
          CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
          QADHF(1,NMST+I)=TEMP4(1,1)
          QADHF(4,NMST+I)=TEMP4(2,1)
          QADHF(2,NMST+I)=TEMP4(3,1)
          QADHF(5,NMST+I)=TEMP4(4,1)
          QADHF(6,NMST+I)=TEMP4(5,1)
          QADHF(3,NMST+I)=TEMP4(6,1)
 602    CONTINUE
C
      END IF
C     }}} Jay
C
C
C ROTATING OCTUPOLE MOMENTS
C
      DO 650 I=1,NPTS
        TEMP2(1) = EFOCT(1,I+NMST)
        TEMP2(2) = EFOCT(4,I+NMST)
        TEMP2(3) = EFOCT(6,I+NMST)
        TEMP2(4) = EFOCT(5,I+NMST)
        TEMP2(5) = EFOCT(10,I+NMST)
        TEMP2(6) = EFOCT(8,I+NMST)
        CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
C
        TEMP2(1)=EFOCT(4,I+NMST)
        TEMP2(2)=EFOCT(6,I+NMST)
        TEMP2(3)=EFOCT(2,I+NMST)
        TEMP2(4)=EFOCT(10,I+NMST)
        TEMP2(5)=EFOCT(7,I+NMST)
        TEMP2(6)=EFOCT(9,I+NMST)
        CALL TFTRI(TEMP4(1,2),TEMP2,T3T,TEMP3,3,3,3)
C
        TEMP2(1)=EFOCT(5,I+NMST)
        TEMP2(2)=EFOCT(10,I+NMST)
        TEMP2(3)=EFOCT(7,I+NMST)
        TEMP2(4)=EFOCT(8,I+NMST)
        TEMP2(5)=EFOCT(9,I+NMST)
        TEMP2(6)=EFOCT(3,I+NMST)
        CALL TFTRI(TEMP4(1,3),TEMP2,T3T,TEMP3,3,3,3)
C
        CALL MRARBR(TEMP4,6,6,3,T3T,3,3,TEMP5,6)
        EFOCT(1,NMST+I)=TEMP5(1,1)
        EFOCT(2,NMST+I)=TEMP5(3,2)
        EFOCT(3,NMST+I)=TEMP5(6,3)
        EFOCT(4,NMST+I)=TEMP5(2,1)
        EFOCT(5,NMST+I)=TEMP5(4,1)
        EFOCT(6,NMST+I)=TEMP5(3,1)
        EFOCT(7,NMST+I)=TEMP5(5,2)
        EFOCT(8,NMST+I)=TEMP5(6,1)
        EFOCT(9,NMST+I)=TEMP5(6,2)
        EFOCT(10,NMST+I)=TEMP5(5,1)
  650 CONTINUE
C
C  ROTATING POLARIZABILITY-COORDINATES AND TENSOR
C
      IPPTS=NPPTS(N)
      DO 670 I=1,IPPTS
        DO 660 J=1,3
          TEMP(J,I) = EFP(J,I+NPST) - EFCENT(J,N)
  660   CONTINUE
  670 CONTINUE
C
      CALL MRARBR(T3,3,3,3,TEMP,3,IPPTS,EFP(1,NPST+1),3)
C
      DO 770 I=1,IPPTS
        DO 760 J=1,3
          EFP(J,I+NPST) = EFP(J,I+NPST) + EFCENT(J,N)
  760   CONTINUE
  770 CONTINUE
C
      IF(PBCWRI) THEN
      DO 679 I = 1,IPPTS
        EFPPB(1,I+NPST) = EFP(1,I+NPST) + XDELPB
        EFPPB(2,I+NPST) = EFP(2,I+NPST) + YDELPB
        EFPPB(3,I+NPST) = EFP(3,I+NPST) + ZDELPB
  679 CONTINUE
      ENDIF
C
      DO 680 I = 1,IPPTS
        EFP(1,I+NPST) = EFP(1,I+NPST) + XDEL
        EFP(2,I+NPST) = EFP(2,I+NPST) + YDEL
        EFP(3,I+NPST) = EFP(3,I+NPST) + ZDEL
  680 CONTINUE
C
      CALL VCLR(TEMP3,1,3)
      DO 690 I=1,IPPTS
        T1(1,1) = EFPOL(1,I+NPST)
        T1(1,2) = EFPOL(4,I+NPST)
        T1(1,3) = EFPOL(5,I+NPST)
        T1(2,1) = EFPOL(7,I+NPST)
        T1(2,2) = EFPOL(2,I+NPST)
        T1(2,3) = EFPOL(6,I+NPST)
        T1(3,1) = EFPOL(8,I+NPST)
        T1(3,2) = EFPOL(9,I+NPST)
        T1(3,3) = EFPOL(3,I+NPST)
        CALL MRARBR(T3,3,3,3,T1,3,3,T2,3)
        CALL VCLR(T1,1,9)
        CALL MRARBR(T2,3,3,3,T3T,3,3,T1,3)
        EFPOL(1,NPST+I) = T1(1,1)
        EFPOL(4,NPST+I) = T1(1,2)
        EFPOL(5,NPST+I) = T1(1,3)
        EFPOL(7,NPST+I) = T1(2,1)
        EFPOL(2,NPST+I) = T1(2,2)
        EFPOL(6,NPST+I) = T1(2,3)
        EFPOL(8,NPST+I) = T1(3,1)
        EFPOL(9,NPST+I) = T1(3,2)
        EFPOL(3,NPST+I) = T1(3,3)
  690 CONTINUE
C
C     Jay {{{
      IF(IEFP2.EQ.1) THEN
        CALL VCLR(TEMP3,1,3)
        DO 691 I=1,IPPTS
          T1(1,1) = EFPOL2(1,I+NPST)
          T1(1,2) = EFPOL2(4,I+NPST)
          T1(1,3) = EFPOL2(5,I+NPST)
          T1(2,1) = EFPOL2(7,I+NPST)
          T1(2,2) = EFPOL2(2,I+NPST)
          T1(2,3) = EFPOL2(6,I+NPST)
          T1(3,1) = EFPOL2(8,I+NPST)
          T1(3,2) = EFPOL2(9,I+NPST)
          T1(3,3) = EFPOL2(3,I+NPST)
          CALL MRARBR(T3,3,3,3,T1,3,3,T2,3)
          CALL VCLR(T1,1,9)
          CALL MRARBR(T2,3,3,3,T3T,3,3,T1,3)
          EFPOL2(1,NPST+I) = T1(1,1)
          EFPOL2(4,NPST+I) = T1(1,2)
          EFPOL2(5,NPST+I) = T1(1,3)
          EFPOL2(7,NPST+I) = T1(2,1)
          EFPOL2(2,NPST+I) = T1(2,2)
          EFPOL2(6,NPST+I) = T1(2,3)
          EFPOL2(8,NPST+I) = T1(3,1)
          EFPOL2(9,NPST+I) = T1(3,2)
          EFPOL2(3,NPST+I) = T1(3,3)
  691   CONTINUE
      END IF
C     }}} Jay
C
C IVANA ADDING DYNAMIC POLARIZABILITY
C  ROTATING DYNAMIC POLARIZABILITY-COORDINATES AND TENSOR
C
      IDPPTS=NDPPTS(N)
      DO 671 I=1,IDPPTS
        DO 661 J=1,3
          DTEMP(J,I) = EFDP(J,I+NDPST) - EFCENT(J,N)
  661   CONTINUE
  671 CONTINUE
C
      CALL MRARBR(DT3,3,3,3,DTEMP,3,IDPPTS,EFDP(1,NDPST+1),3)
C
      DO 771 I=1,IDPPTS
        DO 761 J=1,3
          EFDP(J,I+NDPST) = EFDP(J,I+NDPST) + EFCENT(J,N)
  761   CONTINUE
  771 CONTINUE
C
      IF(PBCWRI) THEN
      DO 678 I = 1,IDPPTS
        EFDPPB(1,I+NDPST) = EFDP(1,I+NDPST) + XDELPB
        EFDPPB(2,I+NDPST) = EFDP(2,I+NDPST) + YDELPB
        EFDPPB(3,I+NDPST) = EFDP(3,I+NDPST) + ZDELPB
  678 CONTINUE
      ENDIF
C
      DO 681 I = 1,IDPPTS
        EFDP(1,I+NDPST) = EFDP(1,I+NDPST) + XDEL
        EFDP(2,I+NDPST) = EFDP(2,I+NDPST) + YDEL
        EFDP(3,I+NDPST) = EFDP(3,I+NDPST) + ZDEL
  681 CONTINUE
C
      CALL VCLR(TEMP3,1,3)
C
      DO 692 I=1,IDPPTS
        T1(1,1) = EFDPOL(1,I+NDPST)
        T1(1,2) = EFDPOL(4,I+NDPST)
        T1(1,3) = EFDPOL(5,I+NDPST)
        T1(2,1) = EFDPOL(7,I+NDPST)
        T1(2,2) = EFDPOL(2,I+NDPST)
        T1(2,3) = EFDPOL(6,I+NDPST)
        T1(3,1) = EFDPOL(8,I+NDPST)
        T1(3,2) = EFDPOL(9,I+NDPST)
        T1(3,3) = EFDPOL(3,I+NDPST)
        CALL MRARBR(T3,3,3,3,T1,3,3,T2,3)
        CALL MRARBR(T2,3,3,3,T3T,3,3,T1,3)
        EFDPOL(1,NDPST+I) = T1(1,1)
        EFDPOL(4,NDPST+I) = T1(1,2)
        EFDPOL(5,NDPST+I) = T1(1,3)
        EFDPOL(7,NDPST+I) = T1(2,1)
        EFDPOL(2,NDPST+I) = T1(2,2)
        EFDPOL(6,NDPST+I) = T1(2,3)
        EFDPOL(8,NDPST+I) = T1(3,1)
        EFDPOL(9,NDPST+I) = T1(3,2)
        EFDPOL(3,NDPST+I) = T1(3,3)
  692 CONTINUE
C
C
      IRPTS = NRPTS(N)
      DO 720 I=1,IRPTS
        DO 710 J=1,3
          TEMP(J,I) = CREP(J,I+NRST) - EFCENT(J,N)
  710   CONTINUE
  720 CONTINUE
C
      CALL MRARBR(T3,3,3,3,TEMP,3,IRPTS,CREP(1,NRST+1),3)
C
      DO 820 I=1,IRPTS
        DO 810 J=1,3
          CREP(J,I+NRST) = CREP(J,I+NRST) + EFCENT(J,N)
  810   CONTINUE
  820 CONTINUE
C
      IF(PBCWRI) THEN
      DO 739 I = 1,IRPTS
        CREPPB(1,I+NRST) = CREP(1,I+NRST) + XDELPB
        CREPPB(2,I+NRST) = CREP(2,I+NRST) + YDELPB
        CREPPB(3,I+NRST) = CREP(3,I+NRST) + ZDELPB
  739 CONTINUE
      ENDIF
C
      DO 740 I = 1,IRPTS
        CREP(1,I+NRST) = CREP(1,I+NRST) + XDEL
        CREP(2,I+NRST) = CREP(2,I+NRST) + YDEL
        CREP(3,I+NRST) = CREP(3,I+NRST) + ZDEL
C        WRITE(6,FMT='(3F9.5)') (CREP(J,I+NRST), J=1,3)
  740 CONTINUE
C
C     Jay {{{
C
      IF(IEFP2.EQ.1) THEN
        IDPTS = NDPTS(N)
        DO 830 I=1,IDPTS
          DO 825 J=1,3
            TEMP(J,I) = CDIS(J,I+NDST) - EFCENT(J,N)
  825     CONTINUE
  830   CONTINUE
C
        CALL MRARBR(T3,3,3,3,TEMP,3,IDPTS,CDIS(1,NDST+1),3)
C
        DO 840 I=1,IDPTS
          DO 835 J=1,3
            CDIS(J,I+NDST) = CDIS(J,I+NDST) + EFCENT(J,N)
  835     CONTINUE
  840   CONTINUE
C
        DO 845 I = 1,IDPTS
          CDIS(1,I+NDST) = CDIS(1,I+NDST) + XDEL
          CDIS(2,I+NDST) = CDIS(2,I+NDST) + YDEL
          CDIS(3,I+NDST) = CDIS(3,I+NDST) + ZDEL
C          WRITE(6,FMT='(3F9.5)') (CDIS(J,I+NDST), J=1,3)
  845   CONTINUE
C
      END IF
C
C     }}} Jay
C
C  PAULI REPULSION INFORMATION
C
      INAT = NATEF(ISET(N))
      DO I=1,INAT
         DO J=1,3
            TEMP(J,I) = PRCORD(J,I+NATST) - EFCENT(J,N)
         END DO
      END DO
      CALL MRARBR(T3,3,3,3,TEMP,3,INAT,PRCORD(1,NATST+1),3)
C
C     if(pbcwri) then
C      do i=1,inat
C         prcopb(1,i+natst) = prcord(1,i+natst) + xdelpb + efcent(1,n)
C         prcopb(2,i+natst) = prcord(2,i+natst) + ydelpb + efcent(2,n)
C         prcopb(3,i+natst) = prcord(3,i+natst) + zdelpb + efcent(3,n)
C      end do
C      endif
C
      DO I=1,INAT
         PRCORD(1,I+NATST) = PRCORD(1,I+NATST) + XDEL + EFCENT(1,N)
         PRCORD(2,I+NATST) = PRCORD(2,I+NATST) + YDEL + EFCENT(2,N)
         PRCORD(3,I+NATST) = PRCORD(3,I+NATST) + ZDEL + EFCENT(3,N)
      END DO
C
      IMO = NORB(N)
      DO NSHL = 1,NSHELL(ISET(N))
         ITYP = KMAX(NSHL,ISET(N)) - KMIN(NSHL,ISET(N)) + 1
         LOC = KLOC(NSHL,ISET(N))
C
C
C  S: DO NOTHING
C
C  P:
C
         IF(ITYP.EQ.3) THEN
            DO I = 1,IMO
               TEMP3(1) = PROVEC(LOC,NORST+I)
               TEMP3(2) = PROVEC(LOC+1,NORST+I)
               TEMP3(3) = PROVEC(LOC+2,NORST+I)
               CALL MRARBR(T3,3,3,3,TEMP3,3,1,TEMP4(1,1),3)
               PROVEC(LOC,NORST+I) = TEMP4(1,1)
               PROVEC(LOC+1,NORST+I) = TEMP4(2,1)
               PROVEC(LOC+2,NORST+I) = TEMP4(3,1)
            END DO
         END IF
C  L:
         IF(ITYP.EQ.4) THEN
            DO I = 1,IMO
               TEMP3(1) = PROVEC(LOC+1,NORST+I)
               TEMP3(2) = PROVEC(LOC+2,NORST+I)
               TEMP3(3) = PROVEC(LOC+3,NORST+I)
               CALL MRARBR(T3,3,3,3,TEMP3,3,1,TEMP4(1,1),3)
               PROVEC(LOC+1,NORST+I) = TEMP4(1,1)
               PROVEC(LOC+2,NORST+I) = TEMP4(2,1)
               PROVEC(LOC+3,NORST+I) = TEMP4(3,1)
            END DO
         END IF
C  D:
C              THE NORMALIZATION FOR XY,XZ,YZ GAUSSIANS IS SQRT(3)
C              LARGER THAN XX,YY,ZZ. FACTOR OF TWO ALLOWS SUMMATION
C              OVER ALL INDICES RATHER THAN I,J<=I
C
         IF(ITYP.EQ.6) THEN
            DO I = 1,IMO
               TEMP2(1) = PROVEC(LOC,NORST+I)
               TEMP2(2) = PROVEC(LOC+3,NORST+I)/CONST
               TEMP2(3) = PROVEC(LOC+1,NORST+I)
               TEMP2(4) = PROVEC(LOC+4,NORST+I)/CONST
               TEMP2(5) = PROVEC(LOC+5,NORST+I)/CONST
               TEMP2(6) = PROVEC(LOC+2,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               PROVEC(LOC,NORST+I) = TEMP4(1,1)
               PROVEC(LOC+3,NORST+I) = TEMP4(2,1)*CONST
               PROVEC(LOC+1,NORST+I) = TEMP4(3,1)
               PROVEC(LOC+4,NORST+I) = TEMP4(4,1)*CONST
               PROVEC(LOC+5,NORST+I) = TEMP4(5,1)*CONST
               PROVEC(LOC+2,NORST+I) = TEMP4(6,1)
            END DO
         END IF
C  F:
         IF(ITYP.EQ.10) THEN
            DO I = 1,IMO
               PROVEC(LOC+3,NORST+I) = PROVEC(LOC+3,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+4,NORST+I) = PROVEC(LOC+4,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+5,NORST+I) = PROVEC(LOC+5,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+6,NORST+I) = PROVEC(LOC+6,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+7,NORST+I) = PROVEC(LOC+7,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+8,NORST+I) = PROVEC(LOC+8,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+9,NORST+I) = PROVEC(LOC+9,NORST+I)*
     *                                 SQRT5/THREE
               PROVEC(LOC+9,NORST+I) = PROVEC(LOC+9,NORST+I)*
     *                                 SQRT3/TWO
C
               TEMP2(1) = PROVEC(LOC  ,NORST+I)
               TEMP2(2) = PROVEC(LOC+3,NORST+I)
               TEMP2(3) = PROVEC(LOC+5,NORST+I)
               TEMP2(4) = PROVEC(LOC+4,NORST+I)
               TEMP2(5) = PROVEC(LOC+9,NORST+I)
               TEMP2(6) = PROVEC(LOC+7,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
C
               TEMP2(1) = PROVEC(LOC+3,NORST+I)
               TEMP2(2) = PROVEC(LOC+5,NORST+I)
               TEMP2(3) = PROVEC(LOC+1,NORST+I)
               TEMP2(4) = PROVEC(LOC+9,NORST+I)
               TEMP2(5) = PROVEC(LOC+6,NORST+I)
               TEMP2(6) = PROVEC(LOC+8,NORST+I)
               CALL TFTRI(TEMP4(1,2),TEMP2,T3T,TEMP3,3,3,3)
C
               TEMP2(1) = PROVEC(LOC+4,NORST+I)
               TEMP2(2) = PROVEC(LOC+9,NORST+I)
               TEMP2(3) = PROVEC(LOC+6,NORST+I)
               TEMP2(4) = PROVEC(LOC+7,NORST+I)
               TEMP2(5) = PROVEC(LOC+8,NORST+I)
               TEMP2(6) = PROVEC(LOC+2,NORST+I)
               CALL TFTRI(TEMP4(1,3),TEMP2,T3T,TEMP3,3,3,3)
C
               CALL MRARBR(TEMP4,6,6,3,T3T,3,3,TEMP5,6)
C
               PROVEC(LOC  ,NORST+I) = TEMP5(1,1)
               PROVEC(LOC+1,NORST+I) = TEMP5(3,2)
               PROVEC(LOC+2,NORST+I) = TEMP5(6,3)
               PROVEC(LOC+3,NORST+I) = TEMP5(2,1)/SQRT5*THREE
               PROVEC(LOC+4,NORST+I) = TEMP5(4,1)/SQRT5*THREE
               PROVEC(LOC+5,NORST+I) = TEMP5(3,1)/SQRT5*THREE
               PROVEC(LOC+6,NORST+I) = TEMP5(5,2)/SQRT5*THREE
               PROVEC(LOC+7,NORST+I) = TEMP5(6,1)/SQRT5*THREE
               PROVEC(LOC+8,NORST+I) = TEMP5(6,2)/SQRT5*THREE
               PROVEC(LOC+9,NORST+I) = TEMP5(5,1)/SQRT5*THREE/
     *                                 SQRT3*TWO
C
            END DO
         END IF
C  G:
C
         IF(ITYP.EQ.15) THEN
            DO I = 1,IMO
               PROVEC(LOC+ 3,NORST+I) = PROVEC(LOC+ 3,NORST+I)*
     *                                  SQRT7/FOUR
               PROVEC(LOC+ 4,NORST+I) = PROVEC(LOC+ 4,NORST+I)*
     *                                  SQRT7/FOUR
               PROVEC(LOC+ 5,NORST+I) = PROVEC(LOC+ 5,NORST+I)*
     *                                  SQRT7/FOUR
               PROVEC(LOC+ 6,NORST+I) = PROVEC(LOC+ 6,NORST+I)*
     *                                  SQRT7/FOUR
               PROVEC(LOC+ 7,NORST+I) = PROVEC(LOC+ 7,NORST+I)*
     *                                  SQRT7/FOUR
               PROVEC(LOC+ 8,NORST+I) = PROVEC(LOC+ 8,NORST+I)*
     *                                  SQRT7/FOUR
               PROVEC(LOC+ 9,NORST+I) = PROVEC(LOC+ 9,NORST+I)*SQRT7*
     *                                  SQRT5/SQRT3/SIX
               PROVEC(LOC+10,NORST+I) = PROVEC(LOC+10,NORST+I)*SQRT7*
     *                                  SQRT5/SQRT3/SIX
               PROVEC(LOC+11,NORST+I) = PROVEC(LOC+11,NORST+I)*SQRT7*
     *                                  SQRT5/SQRT3/SIX
               PROVEC(LOC+12,NORST+I) = PROVEC(LOC+12,NORST+I)*SQRT7*
     *                                  SQRT5/TWELVE
               PROVEC(LOC+13,NORST+I) = PROVEC(LOC+13,NORST+I)*SQRT7*
     *                                  SQRT5/TWELVE
               PROVEC(LOC+14,NORST+I) = PROVEC(LOC+14,NORST+I)*SQRT7*
     *                                  SQRT5/TWELVE
C
C              PERFORM THE FIRST AND SECOND ROTATION
C
               TEMP2(1) = PROVEC(LOC   ,NORST+I)
               TEMP2(2) = PROVEC(LOC+ 3,NORST+I)
               TEMP2(3) = PROVEC(LOC+ 9,NORST+I)
               TEMP2(4) = PROVEC(LOC+ 4,NORST+I)
               TEMP2(5) = PROVEC(LOC+12,NORST+I)
               TEMP2(6) = PROVEC(LOC+10,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXX(1) = TEMP4(1,1)
               TPXX(2) = TEMP4(2,1)
               TPXX(3) = TEMP4(3,1)
               TPXX(4) = TEMP4(4,1)
               TPXX(5) = TEMP4(5,1)
               TPXX(6) = TEMP4(6,1)
C
               TEMP2(1) = PROVEC(LOC+ 3,NORST+I)
               TEMP2(2) = PROVEC(LOC+ 9,NORST+I)
               TEMP2(3) = PROVEC(LOC+ 5,NORST+I)
               TEMP2(4) = PROVEC(LOC+12,NORST+I)
               TEMP2(5) = PROVEC(LOC+13,NORST+I)
               TEMP2(6) = PROVEC(LOC+14,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXY(1) = TEMP4(1,1)
               TPXY(2) = TEMP4(2,1)
               TPXY(3) = TEMP4(3,1)
               TPXY(4) = TEMP4(4,1)
               TPXY(5) = TEMP4(5,1)
               TPXY(6) = TEMP4(6,1)
C
               TEMP2(1) = PROVEC(LOC+ 9,NORST+I)
               TEMP2(2) = PROVEC(LOC+ 5,NORST+I)
               TEMP2(3) = PROVEC(LOC+ 1,NORST+I)
               TEMP2(4) = PROVEC(LOC+13,NORST+I)
               TEMP2(5) = PROVEC(LOC+ 6,NORST+I)
               TEMP2(6) = PROVEC(LOC+11,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYY(1) = TEMP4(1,1)
               TPYY(2) = TEMP4(2,1)
               TPYY(3) = TEMP4(3,1)
               TPYY(4) = TEMP4(4,1)
               TPYY(5) = TEMP4(5,1)
               TPYY(6) = TEMP4(6,1)
C
               TEMP2(1) = PROVEC(LOC+ 4,NORST+I)
               TEMP2(2) = PROVEC(LOC+12,NORST+I)
               TEMP2(3) = PROVEC(LOC+13,NORST+I)
               TEMP2(4) = PROVEC(LOC+10,NORST+I)
               TEMP2(5) = PROVEC(LOC+14,NORST+I)
               TEMP2(6) = PROVEC(LOC+ 7,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXZ(1) = TEMP4(1,1)
               TPXZ(2) = TEMP4(2,1)
               TPXZ(3) = TEMP4(3,1)
               TPXZ(4) = TEMP4(4,1)
               TPXZ(5) = TEMP4(5,1)
               TPXZ(6) = TEMP4(6,1)
C
               TEMP2(1) = PROVEC(LOC+12,NORST+I)
               TEMP2(2) = PROVEC(LOC+13,NORST+I)
               TEMP2(3) = PROVEC(LOC+ 6,NORST+I)
               TEMP2(4) = PROVEC(LOC+14,NORST+I)
               TEMP2(5) = PROVEC(LOC+11,NORST+I)
               TEMP2(6) = PROVEC(LOC+ 8,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYZ(1) = TEMP4(1,1)
               TPYZ(2) = TEMP4(2,1)
               TPYZ(3) = TEMP4(3,1)
               TPYZ(4) = TEMP4(4,1)
               TPYZ(5) = TEMP4(5,1)
               TPYZ(6) = TEMP4(6,1)
C
               TEMP2(1) = PROVEC(LOC+10,NORST+I)
               TEMP2(2) = PROVEC(LOC+14,NORST+I)
               TEMP2(3) = PROVEC(LOC+11,NORST+I)
               TEMP2(4) = PROVEC(LOC+ 7,NORST+I)
               TEMP2(5) = PROVEC(LOC+ 8,NORST+I)
               TEMP2(6) = PROVEC(LOC+ 2,NORST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPZZ(1) = TEMP4(1,1)
               TPZZ(2) = TEMP4(2,1)
               TPZZ(3) = TEMP4(3,1)
               TPZZ(4) = TEMP4(4,1)
               TPZZ(5) = TEMP4(5,1)
               TPZZ(6) = TEMP4(6,1)
C
C              PERFORM THE THIRD AND FOURTH ROTATION
C              PERFORM THE THIRD AND FOURTH ROTATION
C
               TEMP2(1) = TPXX(1)
               TEMP2(2) = TPXY(1)
               TEMP2(3) = TPYY(1)
               TEMP2(4) = TPXZ(1)
               TEMP2(5) = TPYZ(1)
               TEMP2(6) = TPZZ(1)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXX2(1) = TEMP4(1,1)
               TPXX2(2) = TEMP4(2,1)
               TPXX2(3) = TEMP4(3,1)
               TPXX2(4) = TEMP4(4,1)
               TPXX2(5) = TEMP4(5,1)
               TPXX2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(2)
               TEMP2(2) = TPXY(2)
               TEMP2(3) = TPYY(2)
               TEMP2(4) = TPXZ(2)
               TEMP2(5) = TPYZ(2)
               TEMP2(6) = TPZZ(2)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXY2(1) = TEMP4(1,1)
               TPXY2(2) = TEMP4(2,1)
               TPXY2(3) = TEMP4(3,1)
               TPXY2(4) = TEMP4(4,1)
               TPXY2(5) = TEMP4(5,1)
               TPXY2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(3)
               TEMP2(2) = TPXY(3)
               TEMP2(3) = TPYY(3)
               TEMP2(4) = TPXZ(3)
               TEMP2(5) = TPYZ(3)
               TEMP2(6) = TPZZ(3)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYY2(1) = TEMP4(1,1)
               TPYY2(2) = TEMP4(2,1)
               TPYY2(3) = TEMP4(3,1)
               TPYY2(4) = TEMP4(4,1)
               TPYY2(5) = TEMP4(5,1)
               TPYY2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(4)
               TEMP2(2) = TPXY(4)
               TEMP2(3) = TPYY(4)
               TEMP2(4) = TPXZ(4)
               TEMP2(5) = TPYZ(4)
               TEMP2(6) = TPZZ(4)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXZ2(1) = TEMP4(1,1)
               TPXZ2(2) = TEMP4(2,1)
               TPXZ2(3) = TEMP4(3,1)
               TPXZ2(4) = TEMP4(4,1)
               TPXZ2(5) = TEMP4(5,1)
               TPXZ2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(5)
               TEMP2(2) = TPXY(5)
               TEMP2(3) = TPYY(5)
               TEMP2(4) = TPXZ(5)
               TEMP2(5) = TPYZ(5)
               TEMP2(6) = TPZZ(5)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYZ2(1) = TEMP4(1,1)
               TPYZ2(2) = TEMP4(2,1)
               TPYZ2(3) = TEMP4(3,1)
               TPYZ2(4) = TEMP4(4,1)
               TPYZ2(5) = TEMP4(5,1)
               TPYZ2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(6)
               TEMP2(2) = TPXY(6)
               TEMP2(3) = TPYY(6)
               TEMP2(4) = TPXZ(6)
               TEMP2(5) = TPYZ(6)
               TEMP2(6) = TPZZ(6)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPZZ2(1) = TEMP4(1,1)
               TPZZ2(2) = TEMP4(2,1)
               TPZZ2(3) = TEMP4(3,1)
               TPZZ2(4) = TEMP4(4,1)
               TPZZ2(5) = TEMP4(5,1)
               TPZZ2(6) = TEMP4(6,1)
C
               PROVEC(LOC   ,NORST+I) = TPXX2(1)
               PROVEC(LOC+ 1,NORST+I) = TPYY2(3)
               PROVEC(LOC+ 2,NORST+I) = TPZZ2(6)
               PROVEC(LOC+ 3,NORST+I) = TPXY2(1)/SQRT7*FOUR
               PROVEC(LOC+ 4,NORST+I) = TPXZ2(1)/SQRT7*FOUR
               PROVEC(LOC+ 5,NORST+I) = TPYY2(2)/SQRT7*FOUR
               PROVEC(LOC+ 6,NORST+I) = TPYZ2(3)/SQRT7*FOUR
               PROVEC(LOC+ 7,NORST+I) = TPXZ2(6)/SQRT7*FOUR
               PROVEC(LOC+ 8,NORST+I) = TPYZ2(6)/SQRT7*FOUR
               PROVEC(LOC+ 9,NORST+I) = TPYY2(1)/SQRT7/SQRT5*
     *                                  SQRT3*SIX
               PROVEC(LOC+10,NORST+I) = TPZZ2(1)/SQRT7/SQRT5*
     *                                  SQRT3*SIX
               PROVEC(LOC+11,NORST+I) = TPZZ2(3)/SQRT7/SQRT5*
     *                                  SQRT3*SIX
               PROVEC(LOC+12,NORST+I) = TPYZ2(1)/SQRT7/SQRT5*
     *                                  TWELVE
               PROVEC(LOC+13,NORST+I) = TPYZ2(2)/SQRT7/SQRT5*
     *                                  TWELVE
               PROVEC(LOC+14,NORST+I) = TPZZ2(2)/SQRT7/SQRT5*
     *                                  TWELVE
            END DO
         END IF
      END DO
C
C      if(pbcwri) then
C      do i=1,imo
C         centpb(1,i+norst) = centcd(1,i+norst) + xdelpb + efcent(1,n)
C         centpb(2,i+norst) = centcd(2,i+norst) + ydelpb + efcent(2,n)
C         centpb(3,i+norst) = centcd(3,i+norst) + zdelpb + efcent(3,n)
C      end do
C      endif
C
      DO I=1,IMO
         DO J=1,3
            TEMP(J,I) = CENTCD(J,I+NORST) - EFCENT(J,N)
         END DO
      END DO
      CALL MRARBR(T3,3,3,3,TEMP,3,IMO,CENTCD(1,NORST+1),3)
      DO I=1,IMO
         CENTCD(1,I+NORST) = CENTCD(1,I+NORST) + XDEL + EFCENT(1,N)
         CENTCD(2,I+NORST) = CENTCD(2,I+NORST) + YDEL + EFCENT(2,N)
         CENTCD(3,I+NORST) = CENTCD(3,I+NORST) + ZDEL + EFCENT(3,N)
      END DO
C  ROTATE CANONICAL VEC
C
      IAO = NPBF(N)
      DO NSHL = 1,NSHELL(ISET(N))
         ITYP = KMAX(NSHL,ISET(N)) - KMIN(NSHL,ISET(N)) + 1
         LOC = KLOC(NSHL,ISET(N))
C
C  S: DO NOTHING
C
C  P:
C
         IF(ITYP.EQ.3) THEN
            DO I = 1,IAO
               TEMP3(1) = CANVEC(LOC,NAOST+I)
               TEMP3(2) = CANVEC(LOC+1,NAOST+I)
               TEMP3(3) = CANVEC(LOC+2,NAOST+I)
               CALL MRARBR(T3,3,3,3,TEMP3,3,1,TEMP4(1,1),3)
               CANVEC(LOC,NAOST+I) = TEMP4(1,1)
               CANVEC(LOC+1,NAOST+I) = TEMP4(2,1)
               CANVEC(LOC+2,NAOST+I) = TEMP4(3,1)
            END DO
         END IF
C  L:
         IF(ITYP.EQ.4) THEN
            DO I = 1,IAO
               TEMP3(1) = CANVEC(LOC+1,NAOST+I)
               TEMP3(2) = CANVEC(LOC+2,NAOST+I)
               TEMP3(3) = CANVEC(LOC+3,NAOST+I)
               CALL MRARBR(T3,3,3,3,TEMP3,3,1,TEMP4(1,1),3)
               CANVEC(LOC+1,NAOST+I) = TEMP4(1,1)
               CANVEC(LOC+2,NAOST+I) = TEMP4(2,1)
               CANVEC(LOC+3,NAOST+I) = TEMP4(3,1)
            END DO
         END IF
C  D:
C              THE NORMALIZATION FOR XY,XZ,YZ GAUSSIANS IS SQRT(3)
C              LARGER THAN XX,YY,ZZ. FACTOR OF TWO ALLOWS SUMMATION
C              OVER ALL INDICES RATHER THAN I,J<=I
C
         IF(ITYP.EQ.6) THEN
            DO I = 1,IAO
               TEMP2(1) = CANVEC(LOC,NAOST+I)
               TEMP2(2) = CANVEC(LOC+3,NAOST+I)/CONST
               TEMP2(3) = CANVEC(LOC+1,NAOST+I)
               TEMP2(4) = CANVEC(LOC+4,NAOST+I)/CONST
               TEMP2(5) = CANVEC(LOC+5,NAOST+I)/CONST
               TEMP2(6) = CANVEC(LOC+2,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               CANVEC(LOC,NAOST+I) = TEMP4(1,1)
               CANVEC(LOC+3,NAOST+I) = TEMP4(2,1)*CONST
               CANVEC(LOC+1,NAOST+I) = TEMP4(3,1)
               CANVEC(LOC+4,NAOST+I) = TEMP4(4,1)*CONST
               CANVEC(LOC+5,NAOST+I) = TEMP4(5,1)*CONST
               CANVEC(LOC+2,NAOST+I) = TEMP4(6,1)
            END DO
         END IF
C  F:
         IF(ITYP.EQ.10) THEN
            DO I = 1,IAO
               CANVEC(LOC+3,NAOST+I) = CANVEC(LOC+3,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+4,NAOST+I) = CANVEC(LOC+4,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+5,NAOST+I) = CANVEC(LOC+5,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+6,NAOST+I) = CANVEC(LOC+6,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+7,NAOST+I) = CANVEC(LOC+7,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+8,NAOST+I) = CANVEC(LOC+8,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+9,NAOST+I) = CANVEC(LOC+9,NAOST+I)*
     *                                 SQRT5/THREE
               CANVEC(LOC+9,NAOST+I) = CANVEC(LOC+9,NAOST+I)*
     *                                 SQRT3/TWO
C
               TEMP2(1) = CANVEC(LOC  ,NAOST+I)
               TEMP2(2) = CANVEC(LOC+3,NAOST+I)
               TEMP2(3) = CANVEC(LOC+5,NAOST+I)
               TEMP2(4) = CANVEC(LOC+4,NAOST+I)
               TEMP2(5) = CANVEC(LOC+9,NAOST+I)
               TEMP2(6) = CANVEC(LOC+7,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
C
               TEMP2(1) = CANVEC(LOC+3,NAOST+I)
               TEMP2(2) = CANVEC(LOC+5,NAOST+I)
               TEMP2(3) = CANVEC(LOC+1,NAOST+I)
               TEMP2(4) = CANVEC(LOC+9,NAOST+I)
               TEMP2(5) = CANVEC(LOC+6,NAOST+I)
               TEMP2(6) = CANVEC(LOC+8,NAOST+I)
               CALL TFTRI(TEMP4(1,2),TEMP2,T3T,TEMP3,3,3,3)
C
               TEMP2(1) = CANVEC(LOC+4,NAOST+I)
               TEMP2(2) = CANVEC(LOC+9,NAOST+I)
               TEMP2(3) = CANVEC(LOC+6,NAOST+I)
               TEMP2(4) = CANVEC(LOC+7,NAOST+I)
               TEMP2(5) = CANVEC(LOC+8,NAOST+I)
               TEMP2(6) = CANVEC(LOC+2,NAOST+I)
               CALL TFTRI(TEMP4(1,3),TEMP2,T3T,TEMP3,3,3,3)
C
               CALL MRARBR(TEMP4,6,6,3,T3T,3,3,TEMP5,6)
C
               CANVEC(LOC  ,NAOST+I) = TEMP5(1,1)
               CANVEC(LOC+1,NAOST+I) = TEMP5(3,2)
               CANVEC(LOC+2,NAOST+I) = TEMP5(6,3)
               CANVEC(LOC+3,NAOST+I) = TEMP5(2,1)/SQRT5*THREE
               CANVEC(LOC+4,NAOST+I) = TEMP5(4,1)/SQRT5*THREE
               CANVEC(LOC+5,NAOST+I) = TEMP5(3,1)/SQRT5*THREE
               CANVEC(LOC+6,NAOST+I) = TEMP5(5,2)/SQRT5*THREE
               CANVEC(LOC+7,NAOST+I) = TEMP5(6,1)/SQRT5*THREE
               CANVEC(LOC+8,NAOST+I) = TEMP5(6,2)/SQRT5*THREE
               CANVEC(LOC+9,NAOST+I) = TEMP5(5,1)/SQRT5*THREE/
     *                                 SQRT3*TWO
C
            END DO
         END IF
C  G:
C
         IF(ITYP.EQ.15) THEN
            DO I = 1,IAO
               CANVEC(LOC+ 3,NAOST+I) = CANVEC(LOC+ 3,NAOST+I)*
     *                                  SQRT7/FOUR
               CANVEC(LOC+ 4,NAOST+I) = CANVEC(LOC+ 4,NAOST+I)*
     *                                  SQRT7/FOUR
               CANVEC(LOC+ 5,NAOST+I) = CANVEC(LOC+ 5,NAOST+I)*
     *                                  SQRT7/FOUR
               CANVEC(LOC+ 6,NAOST+I) = CANVEC(LOC+ 6,NAOST+I)*
     *                                  SQRT7/FOUR
               CANVEC(LOC+ 7,NAOST+I) = CANVEC(LOC+ 7,NAOST+I)*
     *                                  SQRT7/FOUR
               CANVEC(LOC+ 8,NAOST+I) = CANVEC(LOC+ 8,NAOST+I)*
     *                                  SQRT7/FOUR
               CANVEC(LOC+ 9,NAOST+I) = CANVEC(LOC+ 9,NAOST+I)*SQRT7*
     *                                  SQRT5/SQRT3/SIX
               CANVEC(LOC+10,NAOST+I) = CANVEC(LOC+10,NAOST+I)*SQRT7*
     *                                  SQRT5/SQRT3/SIX
               CANVEC(LOC+11,NAOST+I) = CANVEC(LOC+11,NAOST+I)*SQRT7*
     *                                  SQRT5/SQRT3/SIX
               CANVEC(LOC+12,NAOST+I) = CANVEC(LOC+12,NAOST+I)*SQRT7*
     *                                  SQRT5/TWELVE
               CANVEC(LOC+13,NAOST+I) = CANVEC(LOC+13,NAOST+I)*SQRT7*
     *                                  SQRT5/TWELVE
               CANVEC(LOC+14,NAOST+I) = CANVEC(LOC+14,NAOST+I)*SQRT7*
     *                                  SQRT5/TWELVE
C
C              PERFORM THE FIRST AND SECOND ROTATION
C
               TEMP2(1) = CANVEC(LOC   ,NAOST+I)
               TEMP2(2) = CANVEC(LOC+ 3,NAOST+I)
               TEMP2(3) = CANVEC(LOC+ 9,NAOST+I)
               TEMP2(4) = CANVEC(LOC+ 4,NAOST+I)
               TEMP2(5) = CANVEC(LOC+12,NAOST+I)
               TEMP2(6) = CANVEC(LOC+10,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXX(1) = TEMP4(1,1)
               TPXX(2) = TEMP4(2,1)
               TPXX(3) = TEMP4(3,1)
               TPXX(4) = TEMP4(4,1)
               TPXX(5) = TEMP4(5,1)
               TPXX(6) = TEMP4(6,1)
C
               TEMP2(1) = CANVEC(LOC+ 3,NAOST+I)
               TEMP2(2) = CANVEC(LOC+ 9,NAOST+I)
               TEMP2(3) = CANVEC(LOC+ 5,NAOST+I)
               TEMP2(4) = CANVEC(LOC+12,NAOST+I)
               TEMP2(5) = CANVEC(LOC+13,NAOST+I)
               TEMP2(6) = CANVEC(LOC+14,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXY(1) = TEMP4(1,1)
               TPXY(2) = TEMP4(2,1)
               TPXY(3) = TEMP4(3,1)
               TPXY(4) = TEMP4(4,1)
               TPXY(5) = TEMP4(5,1)
               TPXY(6) = TEMP4(6,1)
               TEMP2(1) = CANVEC(LOC+ 9,NAOST+I)
               TEMP2(2) = CANVEC(LOC+ 5,NAOST+I)
               TEMP2(3) = CANVEC(LOC+ 1,NAOST+I)
               TEMP2(4) = CANVEC(LOC+13,NAOST+I)
               TEMP2(5) = CANVEC(LOC+ 6,NAOST+I)
               TEMP2(6) = CANVEC(LOC+11,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYY(1) = TEMP4(1,1)
               TPYY(2) = TEMP4(2,1)
               TPYY(3) = TEMP4(3,1)
               TPYY(4) = TEMP4(4,1)
               TPYY(5) = TEMP4(5,1)
               TPYY(6) = TEMP4(6,1)
C
               TEMP2(1) = CANVEC(LOC+ 4,NAOST+I)
               TEMP2(2) = CANVEC(LOC+12,NAOST+I)
               TEMP2(3) = CANVEC(LOC+13,NAOST+I)
               TEMP2(4) = CANVEC(LOC+10,NAOST+I)
               TEMP2(5) = CANVEC(LOC+14,NAOST+I)
               TEMP2(6) = CANVEC(LOC+ 7,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXZ(1) = TEMP4(1,1)
               TPXZ(2) = TEMP4(2,1)
               TPXZ(3) = TEMP4(3,1)
               TPXZ(4) = TEMP4(4,1)
               TPXZ(5) = TEMP4(5,1)
               TPXZ(6) = TEMP4(6,1)
C
               TEMP2(1) = CANVEC(LOC+12,NAOST+I)
               TEMP2(2) = CANVEC(LOC+13,NAOST+I)
               TEMP2(3) = CANVEC(LOC+ 6,NAOST+I)
               TEMP2(4) = CANVEC(LOC+14,NAOST+I)
               TEMP2(5) = CANVEC(LOC+11,NAOST+I)
               TEMP2(6) = CANVEC(LOC+ 8,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYZ(1) = TEMP4(1,1)
               TPYZ(2) = TEMP4(2,1)
               TPYZ(3) = TEMP4(3,1)
               TPYZ(4) = TEMP4(4,1)
               TPYZ(5) = TEMP4(5,1)
               TPYZ(6) = TEMP4(6,1)
C
               TEMP2(1) = CANVEC(LOC+10,NAOST+I)
               TEMP2(2) = CANVEC(LOC+14,NAOST+I)
               TEMP2(3) = CANVEC(LOC+11,NAOST+I)
               TEMP2(4) = CANVEC(LOC+ 7,NAOST+I)
               TEMP2(5) = CANVEC(LOC+ 8,NAOST+I)
               TEMP2(6) = CANVEC(LOC+ 2,NAOST+I)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPZZ(1) = TEMP4(1,1)
               TPZZ(2) = TEMP4(2,1)
               TPZZ(3) = TEMP4(3,1)
               TPZZ(4) = TEMP4(4,1)
               TPZZ(5) = TEMP4(5,1)
               TPZZ(6) = TEMP4(6,1)
C
C              PERFORM THE THIRD AND FOURTH ROTATION
C
               TEMP2(1) = TPXX(1)
               TEMP2(2) = TPXY(1)
               TEMP2(3) = TPYY(1)
               TEMP2(4) = TPXZ(1)
               TEMP2(5) = TPYZ(1)
               TEMP2(6) = TPZZ(1)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXX2(1) = TEMP4(1,1)
               TPXX2(2) = TEMP4(2,1)
               TPXX2(3) = TEMP4(3,1)
               TPXX2(4) = TEMP4(4,1)
               TPXX2(5) = TEMP4(5,1)
               TPXX2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(2)
               TEMP2(2) = TPXY(2)
               TEMP2(3) = TPYY(2)
               TEMP2(4) = TPXZ(2)
               TEMP2(5) = TPYZ(2)
               TEMP2(6) = TPZZ(2)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXY2(1) = TEMP4(1,1)
               TPXY2(2) = TEMP4(2,1)
               TPXY2(3) = TEMP4(3,1)
               TPXY2(4) = TEMP4(4,1)
               TPXY2(5) = TEMP4(5,1)
               TPXY2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(3)
               TEMP2(2) = TPXY(3)
               TEMP2(3) = TPYY(3)
               TEMP2(4) = TPXZ(3)
               TEMP2(5) = TPYZ(3)
               TEMP2(6) = TPZZ(3)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYY2(1) = TEMP4(1,1)
               TPYY2(2) = TEMP4(2,1)
               TPYY2(3) = TEMP4(3,1)
               TPYY2(4) = TEMP4(4,1)
               TPYY2(5) = TEMP4(5,1)
               TPYY2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(4)
               TEMP2(2) = TPXY(4)
               TEMP2(3) = TPYY(4)
               TEMP2(4) = TPXZ(4)
               TEMP2(5) = TPYZ(4)
               TEMP2(6) = TPZZ(4)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPXZ2(1) = TEMP4(1,1)
               TPXZ2(2) = TEMP4(2,1)
               TPXZ2(3) = TEMP4(3,1)
               TPXZ2(4) = TEMP4(4,1)
               TPXZ2(5) = TEMP4(5,1)
               TPXZ2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(5)
               TEMP2(2) = TPXY(5)
               TEMP2(3) = TPYY(5)
               TEMP2(4) = TPXZ(5)
               TEMP2(5) = TPYZ(5)
               TEMP2(6) = TPZZ(5)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPYZ2(1) = TEMP4(1,1)
               TPYZ2(2) = TEMP4(2,1)
               TPYZ2(3) = TEMP4(3,1)
               TPYZ2(4) = TEMP4(4,1)
               TPYZ2(5) = TEMP4(5,1)
               TPYZ2(6) = TEMP4(6,1)
C
               TEMP2(1) = TPXX(6)
               TEMP2(2) = TPXY(6)
               TEMP2(3) = TPYY(6)
               TEMP2(4) = TPXZ(6)
               TEMP2(5) = TPYZ(6)
               TEMP2(6) = TPZZ(6)
               CALL TFTRI(TEMP4(1,1),TEMP2,T3T,TEMP3,3,3,3)
               TPZZ2(1) = TEMP4(1,1)
               TPZZ2(2) = TEMP4(2,1)
               TPZZ2(3) = TEMP4(3,1)
               TPZZ2(4) = TEMP4(4,1)
               TPZZ2(5) = TEMP4(5,1)
               TPZZ2(6) = TEMP4(6,1)
C
               CANVEC(LOC   ,NAOST+I) = TPXX2(1)
               CANVEC(LOC+ 1,NAOST+I) = TPYY2(3)
               CANVEC(LOC+ 2,NAOST+I) = TPZZ2(6)
               CANVEC(LOC+ 3,NAOST+I) = TPXY2(1)/SQRT7*FOUR
               CANVEC(LOC+ 4,NAOST+I) = TPXZ2(1)/SQRT7*FOUR
               CANVEC(LOC+ 5,NAOST+I) = TPYY2(2)/SQRT7*FOUR
               CANVEC(LOC+ 6,NAOST+I) = TPYZ2(3)/SQRT7*FOUR
               CANVEC(LOC+ 7,NAOST+I) = TPXZ2(6)/SQRT7*FOUR
               CANVEC(LOC+ 8,NAOST+I) = TPYZ2(6)/SQRT7*FOUR
               CANVEC(LOC+ 9,NAOST+I) = TPYY2(1)/SQRT7/SQRT5*
     *                                  SQRT3*SIX
               CANVEC(LOC+10,NAOST+I) = TPZZ2(1)/SQRT7/SQRT5*
     *                                  SQRT3*SIX
               CANVEC(LOC+11,NAOST+I) = TPZZ2(3)/SQRT7/SQRT5*
     *                                  SQRT3*SIX
               CANVEC(LOC+12,NAOST+I) = TPYZ2(1)/SQRT7/SQRT5*
     *                                  TWELVE
               CANVEC(LOC+13,NAOST+I) = TPYZ2(2)/SQRT7/SQRT5*
     *                                  TWELVE
               CANVEC(LOC+14,NAOST+I) = TPZZ2(2)/SQRT7/SQRT5*
     *                                  TWELVE
            END DO
         END IF
      END DO
C
      CALL EFCM
C
      RETURN
      END
C
C*MODULE MDDSPQ   *DECK MDROTT
      SUBROUTINE MDROTT(QW1,QX1,QY1,QZ1,QW2,QX2,QY2,QZ2,T,INTMD)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      INTEGER INTMD
      DIMENSION T(3,3)
C
C        write(6,*) 'in mdrott:'
C        write(6,*) 'q1 = ',qw1,qx1,qy1,qz1
C        write(6,*) 'q2 = ',qw2,qx2,qy2,qz2
C
       IF(INTMD.EQ.1) THEN
        AXX1 = QW1 ** 2 + QX1 ** 2 - QY1 ** 2 - QZ1 ** 2
        AXY1 = 2.0D+00 * ( QX1 * QY1 + QW1 * QZ1 )
        AXZ1 = 2.0D+00 * ( QX1 * QZ1 - QW1 * QY1 )
        AYX1 = 2.0D+00 * ( QX1 * QY1 - QW1 * QZ1 )
        AYY1 = QW1 ** 2 - QX1 ** 2 + QY1 ** 2 - QZ1 ** 2
        AYZ1 = 2.0D+00 * ( QY1 * QZ1 + QW1 * QX1 )
        AZX1 = 2.0D+00 * ( QX1 * QZ1 + QW1 * QY1 )
        AZY1 = 2.0D+00 * ( QY1 * QZ1 - QW1 * QX1 )
        AZZ1 = QW1 ** 2 - QX1 ** 2 - QY1 ** 2 + QZ1 ** 2
C
        AXX2 = QW2 ** 2 + QX2 ** 2 - QY2 ** 2 - QZ2 ** 2
        AXY2 = 2.0D+00 * ( QX2 * QY2 + QW2 * QZ2 )
        AXZ2 = 2.0D+00 * ( QX2 * QZ2 - QW2 * QY2 )
        AYX2 = 2.0D+00 * ( QX2 * QY2 - QW2 * QZ2 )
        AYY2 = QW2 ** 2 - QX2 ** 2 + QY2 ** 2 - QZ2 ** 2
        AYZ2 = 2.0D+00 * ( QY2 * QZ2 + QW2 * QX2 )
        AZX2 = 2.0D+00 * ( QX2 * QZ2 + QW2 * QY2 )
        AZY2 = 2.0D+00 * ( QY2 * QZ2 - QW2 * QX2 )
        AZZ2 = QW2 ** 2 - QX2 ** 2 - QY2 ** 2 + QZ2 ** 2
       ELSE
         AXX1=-QW1**2+QX1**2-QY1**2+QZ1**2
         AXY1= 2.0D+00*(QY1*QZ1-QW1*QX1)
         AXZ1= 2.0D+00*(QX1*QY1+QW1*QZ1)
C
         AYX1=-2.0D+00*(QW1*QX1+QY1*QZ1)
         AYY1= QW1**2-QX1**2-QY1**2+QZ1**2
         AYZ1= 2.0D+00*(QX1*QZ1-QW1*QY1)
C
         AZX1= 2.0D+00*(QX1*QY1-QW1*QZ1)
         AZY1=-2.0D+00*(QW1*QY1+QX1*QZ1)
         AZZ1=-QW1**2-QX1**2+QY1**2+QZ1**2
C
         AXX2=-QW2**2+QX2**2-QY2**2+QZ2**2
         AXY2= 2.0D+00*(QY2*QZ2-QW2*QX2)
         AXZ2= 2.0D+00*(QX2*QY2+QW2*QZ2)
C
         AYX2=-2.0D+00*(QW2*QX2+QY2*QZ2)
         AYY2= QW2**2-QX2**2-QY2**2+QZ2**2
         AYZ2= 2.0D+00*(QX2*QZ2-QW2*QY2)
C
         AZX2= 2.0D+00*(QX2*QY2-QW2*QZ2)
         AZY2=-2.0D+00*(QW2*QY2+QX2*QZ2)
         AZZ2=-QW2**2-QX2**2+QY2**2+QZ2**2
       ENDIF
C
C
C
C        t should be a2 * a1^T
C
C        t(1,1) = axx2*axx1 + axy2*axy1 + axz2*axz1
C        t(1,2) = axx2*ayx1 + axy2*ayy1 + axz2*ayz1
C        t(1,3) = axx2*azx1 + axy2*azy1 + axz2*azz1
C        t(2,1) = ayx2*axx1 + ayy2*axy1 + ayz2*axz1
C        t(2,2) = ayx2*ayx1 + ayy2*ayy1 + ayz2*ayz1
C        t(2,3) = ayx2*azx1 + ayy2*azy1 + ayz2*azz1
C        t(3,1) = azx2*axx1 + azy2*axy1 + azz2*axz1
C        t(3,2) = azx2*ayx1 + azy2*ayy1 + azz2*ayz1
C        t(3,3) = azx2*azx1 + azy2*azy1 + azz2*azz1
C
C        a2^T * a1
C
        T(1,1) = AXX2*AXX1 + AYX2*AYX1 + AZX2*AZX1
        T(1,2) = AXX2*AXY1 + AYX2*AYY1 + AZX2*AZY1
        T(1,3) = AXX2*AXZ1 + AYX2*AYZ1 + AZX2*AZZ1
        T(2,1) = AXY2*AXX1 + AYY2*AYX1 + AZY2*AZX1
        T(2,2) = AXY2*AXY1 + AYY2*AYY1 + AZY2*AZY1
        T(2,3) = AXY2*AXZ1 + AYY2*AYZ1 + AZY2*AZZ1
        T(3,1) = AXZ2*AXX1 + AYZ2*AYX1 + AZZ2*AZX1
        T(3,2) = AXZ2*AXY1 + AYZ2*AYY1 + AZZ2*AZY1
        T(3,3) = AXZ2*AXZ1 + AYZ2*AYZ1 + AZZ2*AZZ1
C
C
        RETURN
        END
C
C*MODULE MDMISC   *DECK MDRPREF
      SUBROUTINE MDRPREF(ITIME)
C        *********************************************************
C        *** Reports to iw various MD quantities      ***
C        *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXPT=2000, MXFRG=1050, MXATM=2000)
C
      DOUBLE PRECISION MDTYP
      LOGICAL GOPARR, DSKWRK, MASWRK
      LOGICAL         NVT
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /FRGNRG/ CCNRG,CDNRG,CQNRG,CONRG,DDNRG,DQNRG,QQNRG,
     *                ESNRG,REPNRG,POLNRG,POLNRG2,DISNRG,DISNRG8
      COMMON /FUNCT / E, EG(3*MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /MDSIM3/ NVT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C-xxx      common /SYSCET/ SYSMAS,SYSCNT(3)
C
C   *******iw = .log and ip = .dat
C
C
        IF ( MASWRK ) THEN
C
            IF(NVT) WRITE(IP,3999)
            WRITE(IP,4000) ITIME
            WRITE(IP,4012) ESNRG*AU2KCAL
            WRITE(IP,4013) POLNRG*AU2KCAL
            IF(IEFP2.EQ.1) THEN
              WRITE(IP,4030) POLNRG2*AU2KCAL
              WRITE(IP,4031) DISNRG*AU2KCAL
              WRITE(IP,4032) DISNRG8*AU2KCAL
            ENDIF
            WRITE(IP,4014) REPNRG*AU2KCAL
            WRITE(IP,4015) E*AU2KCAL
C-XXX-            write(ip,4016) (syscnt(j)*au2ang,j=1,3)
C-XXX-            write(ip,*) 'syscnt(j),au2ang)',
C-XXX-     *                   syscnt(1),syscnt(2),syscnt(3),au2ang
          IF (ITIME.NE.0) THEN
            EKIN = EKINT + EKINR
            WRITE(IP,4027) EKINT
            WRITE(IP,4028) EKINR
            WRITE(IP,4029) EKIN
            WRITE(IP,4039) (EKIN+E*AU2KCAL)
            WRITE(IP,4050) TEMPMD(EKIN)
          ELSE
            RETURN
          ENDIF
        ENDIF
C
        CALL FLSHBF(IP)
C
        RETURN
C       *** report formats
 3999   FORMAT(/,1X,'NVT RUN!!!')
 4000   FORMAT(/,1X,'*** AT T=',I8)
 4012   FORMAT(  5X,'ELEC ENERGY      = ',F20.8,' KCAL/MOL')
 4013   FORMAT(  5X,'POLZ ENERGY      = ',F20.8,' KCAL/MOL')
 4014   FORMAT(  5X,'REP  ENERGY      = ',F20.8,' KCAL/MOL')
 4015   FORMAT(  5X,'POT  ENERGY      = ',F20.8,' KCAL/MOL')
 4027   FORMAT(/,5X,'TRANS    KE      = ',F20.8,' KCAL/MOL')
 4028   FORMAT(  5X,'ROT      KE      = ',F20.8,' KCAL/MOL')
 4029   FORMAT(  5X,'TOT      KE      = ',F20.8,' KCAL/MOL')
 4030   FORMAT(  5X,'POL2 ENERGY      = ',F20.8,' KCAL/MOL')
 4031   FORMAT(  5X,'DISP ENERGY      = ',F20.8,' KCAL/MOL')
 4032   FORMAT(  5X,'DISP8ENERGY      = ',F20.8,' KCAL/MOL')
 4039   FORMAT(/,5X,'TOT  ENERGY      = ',F20.8,' KCAL/MOL')
 4050   FORMAT(/,5X,'TEMPERATURE      = ',F20.8,' K'       )
        END
C
C*MODULE MDMISC   *DECK MDRPRAB
      SUBROUTINE MDRPRAB(ITIME)
C        *********************************************************
C        *** Reports to iw various MD quantities      ***
C        *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      DOUBLE PRECISION MDTYP
      LOGICAL GOPARR, DSKWRK, MASWRK
      LOGICAL NVT
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /FUNCT / E, EG(3*MXATM)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /MDSIM3/ NVT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C   *******iw = .log and ip = .dat
C
        IF ( MASWRK ) THEN
          IF(NVT) WRITE(IP,3999)
          WRITE(IP,4000) ITIME
          WRITE(IP,*) E*AU2KCAL
          WRITE(IP,4015) E*AU2KCAL
          IF (ITIME.NE.0) THEN
            EKIN = EKINT
            WRITE(IP,4029) EKIN
            WRITE(IP,*) (EKIN+E*AU2KCAL)
            WRITE(IP,4039) (EKIN+E*AU2KCAL)
            WRITE(IP,4050) TEMPMD(EKIN)
          ELSE
            RETURN
          ENDIF
       ENDIF
C
        RETURN
C        *** report formats
 3999   FORMAT(/,1X,'NVT RUN!!!')
 4000        FORMAT(/,1X,'*** AT T=',I8)
 4015   FORMAT(5X,'POT. ENERGY           = ',E14.8,' KCAL/MOL')
 4029   FORMAT(5X,'TOT. K.E.             = ',F14.8,' KCAL/MOL')
 4039   FORMAT(/,5X,'TOT. ENERGY       = ',E14.8,' KCAL/MOL')
 4050        FORMAT(/,5X,'TEMPERATURE           = ',F13.2,' K')
        END
C
C*MODULE MDMISC   *DECK MDWRAB
      SUBROUTINE MDWRAB(ITIME)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXATM=2000, MXAO=8192)
C
      LOGICAL READ
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM2/ READ
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C       *********************************************************
C       *** Writes out MD configuration                       ***
C       *********************************************************
C
C--hmn
      IPIRC = 4
C
      WRITE(IPIRC,2000) NAT
      WRITE(IPIRC,2001) ITIME
C
      DO I = 1, NAT
          WRITE(IPIRC,2010) ANAM(I),C(1,I)*AU2ANG, C(2,I)*AU2ANG,
     *                   C(3,I)*AU2ANG
      ENDDO
C
      IF (ITIME.NE.0.OR.READ) THEN
         VFACT = 1.0D-12*AU2ANG
         WRITE(IW,2001) ITIME
         WRITE(IW,3040)
         DO I = 1, NAT
            WRITE(IW,3050) VX(I)*VFACT, VY(I)*VFACT, VZ(I)*VFACT
         ENDDO
         WRITE(IW,*) '   '
      ENDIF
      RETURN
C
 2000 FORMAT(I8)
 2001 FORMAT('SIMULATION STEP = ',I8)
 2010 FORMAT(1X,A8,F14.10,2X,F14.10,2X,F14.10)
 3040 FORMAT(/,5X,'TRANSLATIONAL VELOCITIES (ANG/PS):')
 3050 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9)
      END
C
C*MODULE MDMISC   *DECK MDWREF
      SUBROUTINE MDWREF(ITIME,INTMD)
C      *********************************************************
C      *** Writes out MD configuration                       ***
C      *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXPT=2000, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      CHARACTER*6 FRGNAM
      CHARACTER*8 PTNAM,FRGNME
C
      LOGICAL READ,MBT,MBR,PBCOUT,PBCWRI,QRAND
C
      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FRAGMT/ XCRD(3,MXFRG),YCRD(3,MXFRG),ZCRD(3,MXFRG),
     *                PTNAM(3,MXFRG),FRGNAM(MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM2/ READ
      COMMON /MDSIM4/ NSTEPS
      COMMON /MDSIM5/ MBT,MBR,QRAND
      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
      COMMON /PBCCO/  FRGCPB(3,MXFGPT),EFPPB(3,MXFGPT),
     *                EFDPPB(3,MXFGPT),CREPPB(3,MXFGPT)
      COMMON /PBC/    PBCOUT,PBCWRI
C--hmn
      IPIRC = 4
C
      WRITE(IPIRC,2000) NFRG*3
      WRITE(IPIRC,2001) ITIME
C
      IAT=0
      KMASS=0
      DO I = 1, NFRG
         WRITE(IPIRC,1010) FRGNAM(I)
         DO III=1,NMPTS(I)
            IAT=IAT+1
            IF(FMASS(III+KMASS).GT.ZERO) THEN
                XF = EFC(1,IAT)*AU2ANG
                YF = EFC(2,IAT)*AU2ANG
                ZF = EFC(3,IAT)*AU2ANG
                WRITE(IPIRC,1020) FRGNME(IAT),XF,YF,ZF
            END IF
         END DO
         KMASS=KMASS+NMPTS(I)
      ENDDO
C
      WRITE(IW,1009)
      IAT=0
      KMASS=0
      DO I = 1, NFRG
          WRITE(IW,1010) FRGNAM(I)
          DO III=1,NMPTS(I)
             IAT=IAT+1
             IF(FMASS(III+KMASS).GT.ZERO) THEN
                XF = EFC(1,IAT)*AU2ANG
                YF = EFC(2,IAT)*AU2ANG
                ZF = EFC(3,IAT)*AU2ANG
                WRITE(IW,1020) FRGNME(IAT),XF,YF,ZF
             END IF
          END DO
          KMASS=KMASS+NMPTS(I)
       ENDDO
C
       IF (ITIME.EQ.NSTEPS) THEN
          IF(PBCWRI) THEN
             WRITE(IPIRC,2002)
             IAT=0
             KMASS=0
             DO I = 1, NFRG
                WRITE(IPIRC,1010) FRGNAM(I)
                DO III=1,NMPTS(I)
                   IAT=IAT+1
                   IF(FMASS(III+KMASS).GT.ZERO) THEN
                      XF = FRGCPB(1,IAT)*AU2ANG
                      YF = FRGCPB(2,IAT)*AU2ANG
                      ZF = FRGCPB(3,IAT)*AU2ANG
                      WRITE(IPIRC,1020) FRGNME(IAT),XF,YF,ZF
                   END IF
                END DO
                KMASS=KMASS+NMPTS(I)
             ENDDO
          ENDIF
C
C            appears to be same as above?
C-XXX-       write(iw,1009)
C-XXX-       IAT=0
C-XXX-       KMASS=0
C-XXX-       do i = 1, nfrg
C-XXX-          write(iw,1010) FRGNAM(i)
C-XXX-          DO III=1,NMPTS(I)
C-XXX-              IAT=IAT+1
C-XXX-              IF(FMASS(III+KMASS).GT.ZERO) THEN
C-XXX-                 XF = EFC(1,IAT)*au2ang
C-XXX-                 YF = EFC(2,IAT)*au2ang
C-XXX-                 ZF = EFC(3,IAT)*au2ang
C-XXX-                 WRITE(IW,1020) FRGNME(IAT),XF,YF,ZF
C-XXX-              END IF
C-XXX-          end do
C-XXX-          KMASS=KMASS+NMPTS(I)
C-XXX-       enddo
C
         IF(PBCWRI) THEN
            WRITE(IW,1011)
            IAT=0
            KMASS=0
            DO I= 1,NFRG
               WRITE(IW,1010) FRGNAM(I)
               DO III=1,NMPTS(I)
                  IAT=IAT+1
                  IF(FMASS(III+KMASS).GT.ZERO) THEN
                     XF =FRGCPB(1,IAT)*AU2ANG
                     YF =FRGCPB(2,IAT)*AU2ANG
                     ZF =FRGCPB(3,IAT)*AU2ANG
                     WRITE(IW,1020) FRGNME(IAT),XF,YF,ZF
                  END IF
               END DO
               KMASS=KMASS+NMPTS(I)
            ENDDO
         ENDIF
      ENDIF
C
      IF (ITIME.NE.0.OR.READ.OR.MBT.OR.MBR) THEN
         VFACT = 1.0D-12
         WRITE(IW,2001) ITIME
         WRITE(IW,3040)
         DO I = 1, NFRG
            WRITE(IW,3050) VX(I)*VFACT, VY(I)*VFACT, VZ(I)*VFACT
         ENDDO
C
         WRITE(IW,3080)
         DO I=1,NFRG
            WRITE(IW,3090) QW(I),QX(I),QY(I),QZ(I)
         ENDDO
C
C--hmn    Distinguish between frog and vverlet...
C
         IF(INTMD.EQ.1) THEN
            WRITE(IW,3060)
            DO I = 1, NFRG
               WRITE(IW,3070) OX(I)*VFACT, OY(I)*VFACT, OZ(I)*VFACT
            ENDDO
            WRITE(IW,*) '   '
         ELSE
            WRITE(IW,4000)
            DO I = 1, NFRG
               WRITE(IW,4010) WX(I)*VFACT, WY(I)*VFACT, WZ(I)*VFACT
            ENDDO
C
            WRITE(IW,4020)
            DO I = 1, NFRG
               WRITE(IW,4030) QW1(I),QX1(I),QY1(I),QZ1(I)
            ENDDO
C
            WRITE(IW,4040)
            DO I = 1, NFRG
               WRITE(IW,4050) QW2(I),QX2(I),QY2(I),QZ2(I)
            ENDDO
            WRITE(IW,*) '   '
         ENDIF
      ELSE
         CONTINUE
      ENDIF
C
      CALL FLSHBF(IW)
      CALL FLSHBF(IPIRC)
      RETURN
C
 1009 FORMAT(11X,'CARTESIAN COORDINATES (ANG)')
 1010 FORMAT(1X,'FRAGNAME=',A6)
 1011 FORMAT(11X,'PBC CARTESIAN COORDINATES (ANG)')
 1020 FORMAT(1X,A8,7X,3F18.13)
 2000 FORMAT(I8)
 2001 FORMAT('SIMULATION STEP = ',I8)
 2002 FORMAT('PBC SIMULATION STEP = ',I8)
 3040 FORMAT(/,5X,'TRANSLATIONAL VELOCITIES (BOHR/PS):')
 3050 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9)
 3060 FORMAT(/,5X,'ANGULAR MOMENTA (RAD/PS):')
 3070 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9)
 3080 FORMAT(/,5X,'QUATERNIONS :')
 3090 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9,2X,E15.9)
 4000 FORMAT(/,5X,'ANGULAR VELOCITY (RAD/PS):')
 4010 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9)
 4020 FORMAT(/,5X,'QUATERNION 1ST DERIV. :')
 4030 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9,2X,E15.9)
 4040 FORMAT(/,5X,'QUATERNION 2ND DERIV. :')
 4050 FORMAT(1X,E15.9,2X,E15.9,2X,E15.9,2X,E15.9)
      END
C
C*MODULE MDMISC   *DECK TEMPMD
      DOUBLE PRECISION FUNCTION TEMPMD(EK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MDTYP
C
      PARAMETER (MXFRG=1050, MXATM=2000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
C
      DATA EFP /8HEFP     /
C
C        tempmd = 2.0d+00 * ek / (6.0d+00*nfrg-3.0d+00)
C     *          / 1.987216d-03
C--pre-IBM
C         tempmd = 2.0d+00 * ek / (6.0d+00*nfrg-6.0d+00)
C     *                 / 1.987216d-03
C
      IF(MDTYP.EQ.EFP) THEN
         TEMPMD = 2.0D+00 * EK / (6.0D+00*NFRG) / 1.987216D-03
      ELSE
         TEMPMD = 2.0D+00 * EK / (6.0D+00*NAT)  / 1.987216D-03
      ENDIF
      RETURN
      END
C
C-XXX-c*module mdmisc   *deck molatm
C-XXX-      subroutine molatm
C-XXX-c        *********************************************************
C-XXX-c        *** Converts quaternion representation to cartesian   ***
C-XXX-c        *********************************************************
C-XXX-      implicit double precision (a-h,o-z)
C-XXX-C
C-XXX-      PARAMETER (MXFRG=1050, MXFGPT=12000,
C-XXX-     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-c
C-XXX-      CHARACTER*8 FRGNME
C-XXX-c
C-XXX-      common /atmdat/ vx(mxfrg), vy(mxfrg), vz(mxfrg),
C-XXX-     *                qw(mxfrg), qx(mxfrg), qy(mxfrg), qz(mxfrg),
C-XXX-     *                ox(mxfrg), oy(mxfrg), oz(mxfrg), qw1(mxfrg),
C-XXX-     *                qx1(mxfrg),qy1(mxfrg),qz1(mxfrg),qw2(mxfrg),
C-XXX-     *                qx2(mxfrg),qy2(mxfrg),qz2(mxfrg)
C-XXX-      common /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),
C-XXX-     *                EFATRM(MXFGPT),
C-XXX-     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),
C-XXX-     *                EFBTRM2(MXFGPT),
C-XXX-     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
C-XXX-     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
C-XXX-      common /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
C-XXX-     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),
C-XXX-     *                FRGMI(6,MXFRG),ATORQ(3,MXFRG)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      common /moddat/ amom(3,MXDFG), totmas(MXDFG), totmai(MXDFG),
C-XXX-     *                drx(MXFGPT), dry(MXFGPT), drz(MXFGPT)
C-XXX-c
C-XXX-      k = 0
C-XXX-      do i = 1, nfrg
C-XXX-          qwi = qw(i)
C-XXX-          qxi = qx(i)
C-XXX-          qyi = qy(i)
C-XXX-          qzi = qz(i)
C-XXX-          axx = qwi ** 2 + qxi ** 2 - qyi ** 2 - qzi ** 2
C-XXX-          axy = 2.0d+00 * ( qxi * qyi + qwi * qzi )
C-XXX-          axz = 2.0d+00 * ( qxi * qzi - qwi * qyi )
C-XXX-          ayx = 2.0d+00 * ( qxi * qyi - qwi * qzi )
C-XXX-          ayy = qwi ** 2 - qxi ** 2 + qyi ** 2 - qzi ** 2
C-XXX-          ayz = 2.0d+00 * ( qyi * qzi + qwi * qxi )
C-XXX-          azx = 2.0d+00 * ( qxi * qzi + qwi * qyi )
C-XXX-          azy = 2.0d+00 * ( qyi * qzi - qwi * qxi )
C-XXX-          azz = qwi ** 2 - qxi ** 2 - qyi ** 2 + qzi ** 2
C-XXX-          do ia = 1, 3
C-XXX-            efc(1,k+ia) = efcent(1,i) + axx * drx(ia)
C-XXX-     *                  + ayx * dry(ia) + azx * drz(ia)
C-XXX-            efc(2,k+ia) = efcent(2,i) + axy * drx(ia)
C-XXX-     *                  + ayy * dry(ia) + azy * drz(ia)
C-XXX-            efc(3,k+ia) = efcent(3,i) + axz * drx(ia)
C-XXX-     *                  + ayz * dry(ia) + azz * drz(ia)
C-XXX-          enddo
C-XXX-          k = k + nmpts(i)
C-XXX-        enddo
C-XXX-        return
C-XXX-        end
C-XXX-C
C-XXX-C*MODULE MDMISC   *DECK MDSCAL
C-XXX-      SUBROUTINE MDSCAL(SCALET,SCALER,INTMD)
C-XXX-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-XXX-C
C-XXX-      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C-XXX-      DOUBLE PRECISION MDTYP
C-XXX-C
C-XXX-      COMMON /ATMDAT/ VX(MXFRG), VY(MXFRG), VZ(MXFRG),
C-XXX-     *                QW(MXFRG), QX(MXFRG), QY(MXFRG), QZ(MXFRG),
C-XXX-     *                OX(MXFRG), OY(MXFRG), OZ(MXFRG), QW1(MXFRG),
C-XXX-     *                QX1(MXFRG),QY1(MXFRG),QZ1(MXFRG),QW2(MXFRG),
C-XXX-     *                QX2(MXFRG),QY2(MXFRG),QZ2(MXFRG)
C-XXX-      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
C-XXX-     *                NPPTS(MXFRG),NPTTPT,IEFP,
C-XXX-     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
C-XXX-     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C-XXX-      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
C-XXX-      COMMON /MDVEL / WX(MXFRG),WY(MXFRG),WZ(MXFRG),TQI(3,MXFRG)
C-XXX-C
C-XXX-      DATA EFP/8HEFP     /
C-XXX-C
C-XXX-C        *************************
C-XXX-C        *** Scales velocities ***
C-XXX-C        *************************
C-XXX-C
C-XXX-      IF(MDTYP.EQ.EFP) THEN
C-XXX-         IF(INTMD.EQ.1) THEN
C-XXX-           DO I = 1, NFRG
C-XXX-            VX(I) = VX(I) * SCALET
C-XXX-            VY(I) = VY(I) * SCALET
C-XXX-            VZ(I) = VZ(I) * SCALET
C-XXX-            OX(I) = OX(I) * SCALER
C-XXX-            OY(I) = OY(I) * SCALER
C-XXX-            OZ(I) = OZ(I) * SCALER
C-XXX-           ENDDO
C-XXX-          ELSE
C-XXX-           DO I = 1, NFRG
C-XXX-            VX(I) = VX(I) * SCALET
C-XXX-            VY(I) = VY(I) * SCALET
C-XXX-            VZ(I) = VZ(I) * SCALET
C-XXX-            WX(I) = WX(I) * SCALER
C-XXX-            WY(I) = WY(I) * SCALER
C-XXX-            WZ(I) = WZ(I) * SCALER
C-XXX-           ENDDO
C-XXX-          ENDIF
C-XXX-        ELSE
C-XXX-          DO I = 1, NFRG
C-XXX-            VX(I) = VX(I) * SCALET
C-XXX-            VY(I) = VY(I) * SCALET
C-XXX-            VZ(I) = VZ(I) * SCALET
C-XXX-          ENDDO
C-XXX-        ENDIF
C-XXX-        RETURN
C-XXX-        END
C
C*MODULE MDMISC   *DECK CALCGR
      SUBROUTINE CALCGR
C        *********************************************************
C        *** calculate goo(r), goh(r), ghh(r) (for one pass)   ***
C        *********************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MDTYP
C
      PARAMETER (MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (NBIN=1000)
C
      DIMENSION INDG(3,3)
C
      CHARACTER*8 FRGNME
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /GOFR  / GOFR(NBIN,3), DELR, NUMG, NUMSAM
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
C
      DATA INDG/1, 2, 2, 2, 3, 3, 2, 3, 3 /
C
C        write(6,*) ' in calcgr'
C        write(6,*) ' delr=',delr
C        write(6,*) ' boxli=',boxli
C---hmn
C        call efcm
C---hmn
        N = 0
         DO I = 1, NFRG
C          comx1=efcent(1,i)
C          comy1=efcent(2,i)
C          comz1 = efcent(3,i)
          M = 0
          DO J = 1, I
            M = M + NMPTS(J)
          ENDDO
          DO J = I+1, NFRG
C            comx2 = efcent(1,J)
C            comy2 = efcent(2,J)
C            comz2 = efcent(3,J)
C            comx = comx1 - comx2
C            comy = comy1 - comy2
C            comz = comz1 - comz2
C
C            pbcx = boxl * anint(comx * boxli)
C            pbcy = boxl * anint(comy * boxli)
C            pbcz = boxl * anint(comz * boxli)
            DO K = 1, 3
              DO L = 1, 3
                GRX = EFC( 1, N+K ) - EFC( 1, M+L )
                GRY = EFC( 2, N+K ) - EFC( 2, M+L )
                GRZ = EFC( 3, N+K ) - EFC( 3, M+L )
                GRX = GRX - BOXL * ANINT( GRX * BOXLI )
                GRY = GRY - BOXL * ANINT( GRY * BOXLI )
                GRZ = GRZ - BOXL * ANINT( GRZ * BOXLI )
C                grx = grx - pbcx
C                gry = gry - pbcy
C                grz = grz - pbcz
                GRL = SQRT(GRX*GRX + GRY*GRY + GRZ*GRZ) * AU2ANG
C                if ( grl .le. 1.0d+00 ) then
C                  write(6,*) ' grl = ',grl
C                  write(6,*) ' i = ',i
C                  write(6,*) ' j = ',j
C                  write(6,*) ' n+k = ',n+k
C                  write(6,*) ' m+l = ',m+l
C                  call abrt
C                endif
                IBIN = INT ( GRL / DELR ) + 1
                IF ( IBIN .LE. NUMG ) THEN
                  INDEX = INDG( L, K )
                  GOFR( IBIN, INDEX ) = GOFR( IBIN, INDEX ) + 1.0D+00
                ENDIF
              ENDDO
            ENDDO
            M = M + NMPTS(J)
          ENDDO
          N = N + NMPTS(I)
        ENDDO
C
        NUMSAM = NUMSAM + 1
C
C        write(6,*) ' leaving calcgr'
C
        RETURN
        END
C
C*MODULE MDMISC   *DECK WRIGR
      SUBROUTINE WRIGR
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MDTYP
C
      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (NBIN=1000)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /CONVMD/ AU2KCAL, EK2KCAL, AU2ANG, AU2SEC, AU2AMU
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /GOFR  / GOFR(NBIN,3), DELR, NUMG, NUMSAM
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MDSIM / DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        **************************
C        *** writes g(r)s to iw ***
C        **************************
C
      IF ( .NOT. MASWRK ) RETURN
C
      IF ( NUMSAM .EQ. 0 ) THEN
          WRITE(IW,3000)
          CALL ABRT
      ENDIF
C
      PI = ACOS(-1.0D+00)
      PREFAC = BOXL*BOXL*BOXL/(2.0D+00*PI/3.0D+00)
      PREFAC = PREFAC / (NFRG*NFRG*NUMSAM)
      PREFAC = PREFAC * AU2ANG * AU2ANG * AU2ANG
C
      DO I = 1, 3
          IF ( I .EQ. 1 ) THEN
            FAC2 = PREFAC
            WRITE(IW,4000)
          ELSE IF ( I .EQ. 2 ) THEN
            FAC2 = 0.25D+00 * PREFAC
            WRITE(IW,4010)
          ELSE
            FAC2 = 0.25D+00 * PREFAC
            WRITE(IW,4020)
          ENDIF
          R = 0.0D+00
          RUPPER = R + DELR
          DO J = 1, NUMG
            DIFF = RUPPER * RUPPER * RUPPER - R * R * R
            WRITE(IW,4030) R+DELR*0.5D+00, GOFR(J,I)*FAC2/DIFF
            R = RUPPER
            RUPPER = RUPPER + DELR
          ENDDO
      ENDDO
C
      CALL FLSHBF(IW)
      RETURN
C
 3000 FORMAT(/,1X,'*** TRYING TO WRITE G(R) WITH NUMSAM=0! ***',/)
 4000 FORMAT(/,5X,'GOO(R)')
 4010 FORMAT(/,5X,'GOH(R)')
 4020 FORMAT(/,5X,'GHH(R)')
 4030 FORMAT(1X,F10.3,5X,E16.10)
      END
