c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     #################################################################
c     ##                                                             ##
c     ##  program dynamic  --  run molecular or stochastic dynamics  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "dynamic" computes a molecular dynamics trajectory in any of
c     several statistical mechanical ensembles with optional periodic
c     boundaries and optional coupling to temperature and pressure baths
c     alternatively a stochastic dynamics trajectory can be generated
c
c
      program dynamic
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'couple.i'
      include 'atmtyp.i'
      include 'bath.i'
      include 'bond.i'
      include 'bound.i'
      include 'inform.i'
      include 'iounit.i'
      include 'charge.i'
      include 'mdstuf.i'
      include 'potent.i'
      include 'solute.i'
      include 'stodyn.i'
      include 'usage.i'
      integer istep,nstep
      integer mode,modstep
      real*8 dt,dtdump
      logical exist,query
      character*120 string
c
c
c     set up the structure and molecular mechanics calculation
c

c     time correlation functions, spectra from FT transformations:
      integer*4 nat0,ia,nw0,iprintft,iw,
     1nw,na1,na2,ix,idumm,nat,nat0ai,NAT1,LP
      parameter (nat0=1000,nw0=1500000,nat0ai=500)
      character*10 s10
      logical lexmd,lwin,lexc,lres,lqftry,lten,ltti,ltei,lexp
      integer*4 iii(4,nat0ai),i,iind,ip,ixp,j,jp,jx,k,kp,kx,l
      real*8 xm(nat0),ym(nat0),zm(nat0),xmm(nat0),ymm(nat0),zmm(nat0),
     2tsux(nw0),tsuy(nw0),tsuz(nw0),tsmx(nw0),tsmy(nw0),tsmz(nw0),
     3tcux(nw0),tcuy(nw0),tcuz(nw0),tcmx(nw0),tcmy(nw0),tcmz(nw0),
     2wsux(nw0),wsuy(nw0),wsuz(nw0),wsmx(nw0),wsmy(nw0),wsmz(nw0),
     3wcux(nw0),wcuy(nw0),wcuz(nw0),wcmx(nw0),wcmy(nw0),wcmz(nw0),
     5spec(nw0),stem(nw0),t,w,dbb,
     1dw,wmax,tp,sol,wmin,res,sc,fac,FBOHR,dtau,bohr,pi
      real*8 fcar,xeq,pchgft(nat0),wfmax,wsup,e(3*nat0ai),pol(nat0),
     1fnew(3*nat0ai,3*nat0ai),SMAT(3*nat0ai,3*nat0ai),Q(3*nat0ai),
     2zmass(3*nat0ai),zim(3*nat0ai),tem(3*nat0ai,3*nat0ai),wwin,
     3dtwin,P(nat0ai,3,3),A(nat0ai,3,3),vx,vy,vz,ux,
     4ALPHAI(3*nat0ai,3,3),AAI(3*nat0ai,3,3,3),GI(3*nat0ai,3,3),
     4ALPHA(3*nat0,3,3),   AA(3*nat0,3,3,3),   G(3*nat0,3,3),
     5APTI(nat0ai,3,3),AI(nat0ai,3,3),u(3,3,nat0ai),Xwork(3,nat0),
     2wsat(3,3,nw0),wcat(3,3,nw0),wsgt(3,3,nw0),wcgt(3,3,nw0),
     2wsaat(3,3,3,nw0),wcaat(3,3,3,nw0),sa1,sal0,sal1,sag0,sag1,
     1s180,d180,ydy,ydx,tf,at(3,3),gt(3,3),aat(3,3,3),s1,s2,
     1at1(3,3),at2(3,3),pt1(3,3),pt2(3,3),ut(3,3),rki(3),rd,
     1gs1(3,3,3),as1(3,3,3),au1(3,3,3,3),au2(3,3,3,3),as2(3,3,3),
     1gs2(3,3,3)
      logical laic,lfixed,lfilterdiag
      common /aic/fcar(3*nat0ai,3*nat0ai),xeq(3*nat0ai),laic
c
c     spectral gradient:
      logical lexsp
      integer*4 icurrent
      real*8 gradsp,xold,yold,zold,ensp,spalpha,sumspold,rshold
      common /spegrad/gradsp(3*maxatm),xold(maxatm),yold(maxatm),
     1zold(maxatm),rshold(maxatm),ensp,lexsp,spalpha,sumspold

      call initial
      call getxyz
      call mechanic
c
c     initialize the temperature, pressure and coupling baths
c
      kelvin = 0.0d0
      atmsph = 0.0d0
      isothermal = .false.
      isobaric = .false.
c
c     initialize the simulation length as number of time steps
c
      query = .true.
      call nextarg (string,exist)
      if (exist) then
         read (string,*,err=10,end=10)  nstep
         query = .false.
      end if
   10 continue
      if (query) then
         write (iout,20)
   20    format (/,' Enter the Number of Dynamics Steps to be',
     &              ' Taken :  ',$)
         read (input,30)  nstep
   30    format (i10)
      end if
c
c     get the length of the dynamics time step in picoseconds
c
      dt = -1.0d0
      call nextarg (string,exist)
      if (exist)  read (string,*,err=40,end=40)  dt
   40 continue
      if (dt .le. 0.0d0) then
         write (iout,50)
   50    format (/,' Enter the Time Step Length in Femtoseconds',
     &              ' [1.0] :  ',$)
         read (input,60)  dt
   60    format (f20.0)
      end if
      if (dt .le. 0.0d0)  dt = 1.0d0
      dt = 1.0d-3 * dt
c
c     set bounds on the Berendsen thermostat coupling parameters
c
      tautemp = max(tautemp,dt)
      taupres = max(taupres,dt)
c
c     set the time between trajectory snapshot coordinate dumps
c
      dtdump = -1.0d0
      call nextarg (string,exist)
      if (exist)  read (string,*,err=70,end=70)  dtdump
   70 continue
      if (dtdump .le. 0.0d0) then
         write (iout,80)
   80    format (/,' Enter Time between Dumps in Picoseconds',
     &              ' [0.1] :  ',$)
         read (input,90)  dtdump
   90    format (f20.0)
      end if
      if (dtdump .le. 0.0d0)  dtdump = 0.1d0
      iwrite = nint(dtdump/dt)
c
c     get choice of statistical ensemble for periodic system
c
      if (use_bounds) then
         mode = 0
         call nextarg (string,exist)
         if (exist)  read (string,*,err=100,end=100)  mode
  100    continue
         if (mode.lt.1 .or. mode.gt.4) then
            write (iout,110)
  110       format (/,' Possible Statistical Mechanical Ensembles :',
     &              //,4x,'(1) Microcanonical (NVE)',
     &              /,4x,'(2) Canonical (NVT)',
     &              /,4x,'(3) Isoenthalpic-Isobaric (NPH)',
     &              /,4x,'(4) Isothermal-Isobaric (NPT)',
     &              //,' Enter the Number of the Desired Choice',
     &                 ' [1] :  ',$)
            read (input,120)  mode
  120       format (i10)
         end if
         if (mode.eq.2 .or. mode.eq.4) then
            isothermal = .true.
            kelvin = -1.0d0
            call nextarg (string,exist)
            if (exist)  read (string,*,err=130,end=130)  kelvin
  130       continue
            if (kelvin .le. 0.0d0) then
               write (iout,140)
  140          format (/,' Enter the Desired Temperature in Degrees',
     &                    ' K [298] :  ',$)
               read (input,150)  kelvin
  150          format (f20.0)
            end if
            if (kelvin .le. 0.0d0)  kelvin = 298.0d0
            kelvin0 = kelvin
         end if
         if (mode.eq.3 .or. mode.eq.4) then
            isobaric = .true.
            atmsph = -1.0d0
            call nextarg (string,exist)
            if (exist)  read (string,*,err=160,end=160)  atmsph
  160       continue
            if (atmsph .le. 0.0d0) then
               write (iout,170)
  170          format (/,' Enter the Desired Pressure in Atm',
     &                    ' [1.0] :  ',$)
               read (input,180)  atmsph
  180          format (f20.0)
            end if
            if (atmsph .le. 0.0d0)  atmsph = 1.0d0
         end if
      end if
c
c     get desired temperature for nonperiodic system
c
      if (.not. use_bounds) then
         isothermal = .true.
         kelvin = -1.0d0
         call nextarg (string,exist)
         if (exist)  read (string,*,err=190,end=190)  kelvin
  190    continue
         if (kelvin .le. 0.0d0) then
            write (iout,200)
  200       format (/,' Enter the Desired Temperature in Degrees',
     &                 ' K [298] :  ',$)
            read (input,210)  kelvin
  210       format (f20.0)
         end if
         if (kelvin .le. 0.0d0)  kelvin = 298.0d0
      end if
c
c     initialize any rattle constraints and setup dynamics
c
      call shakeup
      call mdinit
c
c     print out a header line for the dynamics computation
c
      if (integrate .eq. 'VERLET') then
         write (iout,220)
  220    format (/,' Molecular Dynamics Trajectory via',
     &              ' Velocity Verlet Algorithm')
      else if (integrate .eq. 'STOCHASTIC') then
         write (iout,230)
  230    format (/,' Stochastic Dynamics Trajectory via',
     &              ' Velocity Verlet Algorithm')
      else if (integrate .eq. 'RIGIDBODY') then
         write (iout,240)
  240    format (/,' Molecular Dynamics Trajectory via',
     &              ' Rigid Body Algorithm')
      else
         write (iout,250)
  250    format (/,' Molecular Dynamics Trajectory via',
     &              ' Modified Beeman Algorithm')
      end if
c
c     integrate equations of motion to take a time step
c
c     Petr Bour's change: enable correlation calculation here:
c     ############################################################
      bohr=0.529177d0
      pi=4.0d0*datan(1.0d0)
c     au^2 -> debye^2
      dbb=2.541765d0**2
      write(6,*)pi
      inquire(file='COR.OPT',exist=lexmd)
      if(lexmd)then
      lres=.false.
      laic=.false.
      lqftry=.true.
      lfilterdiag=.false.
      tp=2.0d0*pi
      sol=3.0d10
c     dw,wmin,wmax step and frequency limits in cm-1
      if(n.gt.nat0)then
       write(6,*)'too many atoms for FT'
       stop
      endif
      na1=1
      na2=n
      iprintft=1
      wmin=0.0d0
      res=-9.0d0
      wmax=4000.0d0
      dw=1.0d0
      lwin=.false.
      lten=.false.
      ltti=.false.
      ltei=.false.
      wwin=1.0d0
      open(7,file='COR.OPT')
111   read(7,71)s10
71    format(a10)
      write(6,71)s10
      if(s10(1:2).eq.'DW')read(7,*)dw
      if(s10(1:2).eq.'N1')read(7,*)na1,na2
      if(s10(1:2).eq.'RE')read(7,*)res
      if(s10(1:4).eq.'WMIN')read(7,*)wmin
      if(s10(1:4).eq.'WMAX')read(7,*)wmax
      if(s10(1:4).eq.'LQFT')read(7,*)lqftry
      if(s10(1:4).eq.'LAIC')read(7,*)laic
      if(s10(1:4).eq.'TENS')read(7,*)lten
      if(s10(1:4).eq.'TTTI')read(7,*)ltti
      if(s10(1:4).eq.'TENI')read(7,*)ltei
      if(s10(1:2).eq.'IP')read(7,*)iprintft
      if(s10(1:4).eq.'CONT')read(7,*)lres
      if(s10(1:5).eq.'FIXED')read(7,*)lfixed,wfmax,wsup
      if(s10(1:5).eq.'FILTE')read(7,*)lfilterdiag
      if(s10(1:3).eq.'WIN')read(7,*)lwin
      if(s10(1:3).eq.'WWI')read(7,*)wwin
      if(s10(1:3).eq.'END')goto 222
      goto 111
222   close(7)
      if(iprintft.ge.1)write(6,*)'FT options have been set'

      if(laic)then
       write(6,*)'ab initio force field FILE.FC'
       write(6,*)'and equilibrium geometry FILE.X used'
       open(88,file='FILE.X')
       read(88,*)
       read(88,*)nat
       if(nat.ne.n)then
        write(6,*)'number of atoms does not agree'
        stop
       endif
       if(nat.gt.nat0ai)then
        write(6,*)'too many atoms for ab initio FF'
        stop
       endif
       do ia=1,nat
        read(88,*)idumm,(xeq(3*(ia-1)+ix),ix=1,3)
       enddo
       close(88)
       call READFF(nat0ai*3,3*nat,FCAR)
c
c      if desired, replace normal modes with frequency smaller than
c      wfmax by wsup:
       if(lfixed)call projectff(nat0ai*3,3*nat,FCAR,xeq,wfmax,wsup,
     1 mass,e,fnew,SMAT,Q,zmass,ZIM,TEM)
      endif

      do 22 ia=1,n
22    pchgft(ia)=pchg(ia)

      inquire(file='CHARGE.LST',exist=lexc)
      if(lexc)then
       write(6,*)'CHARGE.LST found, partial charges read from it'
       open(99,file='CHARGE.LST')
       read(99,*)
       do 19 ia=1,n
19     read(99,*)pchgft(ia)
       close(99)
      endif

      inquire(file='POLAR.LST',exist=lexp)
      if(lexp)then
       write(6,*)'POLAR.LST found, atomic polarizabilities read from it'
       open(99,file='POLAR.LST')
       read(99,*)
       do 18 ia=1,n
18     read(99,*)pol(ia)
       close(99)
      endif

      if(lten)call READTEN(nat0ai,P,A)
      if(ltti)then
       open(15,FILE='INT.TTT',status='old')
       CALL READTTT(nat0ai,ALPHAI,AAI,GI,NAT1)
       call readatms(15,NAT1,iii,nat0ai)
       close(15)
       if(NAT1.ne.n)call report('Strange number of atoms in INT.TTT')
      endif
      if(ltei)then
       OPEN(15,FILE='INT.TEN',status='old')
       CALL READTENI(nat0ai,APTI,AI,N)
       call readatms(15,N,iii,nat0ai)
       CLOSE(15)
      endif

      inquire(file='FTRY.INP',exist=lexc)
      if(lexc)then
       write(6,*)'FTRY.INP found, masses charges read from it'
       write(6,*)'old masses:'
       write(6,6009)(mass(ia),ia=1,n)
6009   format(6f12.6)
       open(99,file='FTRY.INP')
       read(99,*)
       read(99,*)
       if(lqftry)then
        read(99,6009)(pchgft(ia),ia=1,n)
        write(6,*)'partial charges read from FTRY.INP'
       else
c       dummy mass reading
        read(99,6009)(mass(ia),ia=1,n)
       endif
       read(99,*)
       read(99,6009)(mass(ia),ia=1,n)
       write(6,*)'new masses:'
       write(6,6009)(mass(ia),ia=1,n)
       close(99)
      endif

      dw=dw*tp*sol
      wmax=wmax*tp*sol
      wmin=wmin*tp*sol
      nw=int((wmax-wmin)/dw)
      if(nw.gt.nw0)then
       write(6,*)'too many spectral points'
       stop
      endif
      sc=1.0d0
      dtwin=dt*wwin*3.0d10*1.0d-12
      if(lwin)sc=exp(-dtwin)

      if(iprintft.ge.1)then
       write(6,*)'atoms ',na1,' - ',na2
       write(6,*)nw,' frequencies'
       if(lwin)write(6,*)'exponentially damped window applied'
       write(6,6700)sc
6700   format(' Scaling factor ',f9.6)
       dtau=dt*1.0d-12/2.418D-17
       write(6,6701)dt,dtau
6701   format(f10.5,' ps step = ',f10.4,' au')
c      time step in atomic units:
       write(6,6702)res
6702   format(f10.5,' cm-1 resolution via convolution')
      endif
c
c     initiate frequency coefficients
c     zero out first to avoid compiler messages:
      call ifc(nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux,
     1wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux,
     3tsuy,tcuy,tsuz,tcuz,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0)
      if(lres)then
       call rwfc(1,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux,
     1 wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux,
     3 tsuy,tcuy,tsuz,tcuz,wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0)
       write(6,*)'Restart from disk'
      else
       write(6,*)'Coefficients zeroed-out'
      endif
      do 24 ia=1,n
      xm(ia)=x(ia)
      ym(ia)=y(ia)
      zm(ia)=z(ia)
      xmm(ia)=x(ia)
      ymm(ia)=y(ia)
24    zmm(ia)=z(ia)
      endif
c     ########################################################
      if(lfilterdiag)open(66,file='UXT.PRN')


c
c     Petr Bour's change: calculate spectrum and S-gradient
c     ########################################################
      inquire(file='SPEC.OPT',exist=lexsp)
      if(lexsp)then
       write(6,*)'SPEC.OPT found'
       call cctn_tinker1(n,atomic,x,y,z,n12,i12,maxval,maxatm)
       call loadnewaoptions(x,y,z,atomic,spalpha)
       call ioparspecomp
       open(67,file='DIFF.PRN')
       do 25 JP=1,3*n
       rshold(JP)=0.0d0
25     gradsp(JP)=0.0d0
       sumspold=0.0d0
      endif
c     ########################################################


      do istep = 1, nstep
c
c      Petr Bour's change: calculate spectrum and S-gradient
c      ########################################################
       if(lexsp)then
        icurrent=istep
        call dospectrum(x,y,z,n12,i12,maxatm,maxval,icurrent)
       endif
c      ########################################################
         if (integrate .eq. 'VERLET') then
            call verlet (istep,dt)
         else if (integrate .eq. 'STOCHASTIC') then
            call sdstep (istep,dt)
         else if (integrate .eq. 'RIGIDBODY') then
            call rgdstep (istep,dt)
         else
            call beeman (istep,dt)
         end if
c
c     remove center of mass translation and rotation if needed
c
         modstep = mod(istep,iprint)
         if (modstep.eq.0 .and. nuse.eq.n) then
            if (integrate.ne.'STOCHASTIC' .and.
     &          thermostat.ne.'ANDERSEN')  call mdrest
         end if
c
c     Petr Bour's change: enable correlation calculation here:
c     ########################################################
      if(lexmd)then
c     update the correlation functions:
c
c     current time [s]
      t=dt*dble(istep)*1.0d-12
c     write(6,*)'update ',int(t*1.0d15),' fs' 
c
      if(istep.gt.2)then
c
c      internal coordinate tensors, define the transformation
c      matrices:
       if(ltei.or.ltti)call dtm(iii,n,u,nat0ai,x,y,z)
       
c      get the correlation function at this time
c      A) selected atoms
       if(ltei)then
        DO 14 L=1,NAT
        DO 1411 IP=1,3
        DO 1411 J=1,3
1411    ut(IP,J)=u(IP,J,L)
        DO 141 IP=1,3
        DO 141 JP=1,3
        at1(IP,JP)=AI(L,IP,JP)
141     pt1(IP,JP)=APTI(L,IP,JP)
        DO 142 I=1,3
        DO 142 JP=1,3
	at2(I,JP)=0.0d0
	pt2(I,JP)=0.0d0
	do 142 IP=1,3
	at2(I,JP)=at2(I,JP)+ut(IP,I)*at1(IP,JP)
142     pt2(I,JP)=pt2(I,JP)+ut(IP,I)*pt1(IP,JP)
        DO 14 I=1,3
        DO 14 J=1,3
	s1=0.0d0
	s2=0.0d0
	do 143 JP=1,3
	s1=s1+ut(JP,J)*at2(I,JP)
143     s2=s2+ut(JP,J)*pt2(I,JP)
        A(L,I,J)=s1
14      P(L,I,J)=s2
c
c       Transform to common origin:
        FBOHR=0.25d0/bohr
        DO 2201 L=1,NAT
        DO 2201 J=1,3
        A(L,J,1)=A(L,J,1)+(y(L)*P(L,J,3)-z(L)*P(L,J,2))*FBOHR
        A(L,J,2)=A(L,J,2)+(z(L)*P(L,J,1)-x(L)*P(L,J,3))*FBOHR
2201    A(L,J,3)=A(L,J,3)+(x(L)*P(L,J,2)-y(L)*P(L,J,1))*FBOHR

       endif

       call updt(na1,na2,wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,
     1 wcmz,wsux,wcux,wsuy,wcuy,wsuz,wcuz,t,sc,xm,ym,zm,xmm,ymm,zmm,
     2 bohr,dtau,x,y,z,A,P,nat0ai,ltei,lten,pchgft,ux)
c
c      polarizability from atomic polarizabilities
       if(lexp)then 
        DO 172 L=1,n
        DO 172 IX=1,3
        IIND=3*(L-1)+IX
        DO 172 I=1,3
        DO 172 J=1,3
        DO 177 K=1,3
177     AA(IIND,I,J,K)=0.0d0
        G(IIND,I,J)=0.0d0
172     ALPHA(IIND,I,J)=0.0d0
        DO 17 L=1,n
        DO 17 LP=1,n
	rki(1)=x(LP)-x(L)
	rki(2)=y(LP)-y(L)
	rki(3)=z(LP)-z(L)
	if(LP.ne.L)then
	 rd=dsqrt(rki(1)**2+rki(2)**2+rki(3)**2+1.0d-9)
         DO 173 IX=1,3
         IIND=3*(L-1)+IX
         DO 173 I=1,3
         DO 173 J=1,3
	 s1=5.0d0*rki(I)*rki(J)*rki(IX)
	 if(I.eq.IX)s1=s1-rd*rd*rki(J)
	 if(J.eq.IX)s1=s1-rd*rd*rki(I)
	 if(I.eq.J)s1=s1-rd*rd*rki(IX)
173      ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+6.0d0*pol(L)*pol(LP)*s1/rd**7
        endif
17      continue
        CALL TTTT(3*nat0,ALPHA,AA,G,NAT,2,x,y,z,Xwork,nat0)
       endif

       if(ltti)then
        DO 15 L=1,n
c       rewrite
        DO 1412 IP=1,3
        DO 1412 J=1,3
1412    ut(IP,J)=u(IP,J,L)
        DO 151 IXP=1,3
        IIND=3*(L-1)+IXP
        DO 151 IP=1,3
        DO 151 JP=1,3
        as1(IXP,IP,JP)=ALPHAI(IIND,IP,JP)
151     gs1(IXP,IP,JP)=G(IIND,IP,JP)
c       first index trafo
        DO 154 IX=1,3
        DO 154 IP=1,3
        DO 154 JP=1,3
        as2(IX,IP,JP)=0.0d0
        gs2(IX,IP,JP)=0.0d0
        DO 154 IXP=1,3
        as2(IX,IP,JP)=as2(IX,IP,JP)+ut(IXP,IX)*as1(IXP,IP,JP)
154     gs2(IX,IP,JP)=gs2(IX,IP,JP)+ut(IXP,IX)*gs1(IXP,IP,JP)
c       second index trafo
        DO 152 IX=1,3
        DO 152 I=1,3
        DO 152 JP=1,3
        as1(IX,I,JP)=0.0d0
        gs1(IX,I,JP)=0.0d0
        DO 152 IP=1,3
        as1(IX,I,JP)=as1(IX,I,JP)+ut(IP,I)*as2(IX,IP,JP)
152     gs1(IX,I,JP)=gs1(IX,I,JP)+ut(IP,I)*gs2(IX,IP,JP)
c       last index trafo
        DO 153 IX=1,3
        IIND=3*(L-1)+IX
        DO 153 I=1,3
        DO 153 J=1,3
        ALPHA(IIND,I,J)=0.0d0
        G(IIND,I,J)=0.0d0
        DO 153 JP=1,3
        ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+ut(JP,J)*as1(IX,I,JP)
153     G(IIND,I,J)    =G(IIND,I,J)    +ut(JP,J)*gs1(IX,I,JP)
c       A-rewrite
        DO 16 IXP=1,3
        IIND=3*(L-1)+IXP
        DO 16 IP=1,3
        DO 16 JP=1,3
        DO 16 KP=1,3
16      au1(IXP,IP,JP,KP)=AAI(IIND,IP,JP,KP)
c       first
        DO 161 IX=1,3
        DO 161 IP=1,3
        DO 161 JP=1,3
        DO 161 KP=1,3
        au2(IX,IP,JP,KP)=0.0d0
        DO 161 IXP=1,3
161     au2(IX,IP,JP,KP)=au2(IX,IP,JP,KP)+ut(IXP,IX)*au1(IXP,IP,JP,KP)
c       second
        DO 162 IX=1,3
        DO 162 I =1,3
        DO 162 JP=1,3
        DO 162 KP=1,3
        au1(IX,I,JP,KP)=0.0d0
        DO 162 IP=1,3
162     au1(IX,I,JP,KP)=au1(IX,I,JP,KP)+ut(IP,I)*au2(IX,IP,JP,KP)
c       third
        DO 163 IX=1,3
        DO 163 I =1,3
        DO 163 J =1,3
        DO 163 KP=1,3
        au2(IX,I,J,KP)=0.0d0
        DO 163 JP=1,3
163     au2(IX,I,J,KP)=au2(IX,I,J,KP)+ut(JP,J)*au1(IX,I,JP,KP)
c       fourth
        DO 15 IX=1,3
        IIND=3*(L-1)+IX
        DO 15 I =1,3
        DO 15 J =1,3
        DO 15 K =1,3
        AA(IIND,I,J,K)=0.0d0
        DO 15 KP=1,3
15      AA(IIND,I,J,K)=AA(IIND,I,J,K)+ut(KP,K)*au2(IX,I,J,KP)
        CALL TTTT(3*nat0ai,ALPHA,AA,G,NAT,2,x,y,z,Xwork,nat0)
       endif

c      time polarizability derivatives at this time
       do 2001 ix=1,3
       do 2001 jx=1,3
       at(ix,jx)=0.0d0
       gt(ix,jx)=0.0d0
       do 2001 kx=1,3
2001   aat(ix,jx,kx)=0.0d0
       if(ltti.or.lexp)then
        fac=1.0d0/dtau/bohr
        do 2002 ia=na1,na2
	vx=(x(ia)-xm(ia))*fac
	vy=(y(ia)-ym(ia))*fac
	vz=(z(ia)-zm(ia))*fac
        do 2002 ix=1,3
        do 2002 jx=1,3
        gt(ix,jx)=gt(ix,jx)        +G(3*(ia-1)+1,ix,jx)*vx
     1  +    G(3*(ia-1)+2,ix,jx)*vy+G(3*(ia-1)+3,ix,jx)*vz
        at(ix,jx)=at(ix,jx)        +ALPHA(3*(ia-1)+1,ix,jx)*vx
     1  +ALPHA(3*(ia-1)+2,ix,jx)*vy+ALPHA(3*(ia-1)+3,ix,jx)*vz
        do 2002 kx=1,3
2002    aat(ix,jx,kx)=aat(ix,jx,kx) +AA(3*(ia-1)+1,ix,jx,kx)*vx
     1  +AA(3*(ia-1)+2,ix,jx,kx)*vy+AA(3*(ia-1)+3,ix,jx,kx)*vz
        call update3(wmin,dw,nw,wsat,wcat,wsgt,wcgt,
     1  wsaat,wcaat,at,gt,aat,t,sc,3*nat0ai)
       endif

       if(lfilterdiag)write(66,661)t*1.0d15,ux/1.0d15
  661  format(2f18.6)
c
c      B) all atoms
       call updt(1,n,wmin,dw,nw,tsmx,tcmx,tsmy,tcmy,tsmz,
     1 tcmz,tsux,tcux,tsuy,tcuy,tsuz,tcuz,t,sc,xm,ym,zm,xmm,ymm,zmm,
     2 bohr,dtau,x,y,z,A,P,nat0ai,ltei,lten,pchgft,ux)

      endif
c     istep>2

c     remember last point
      do 1 ia=1,n
      xmm(ia)=xm(ia)
      ymm(ia)=ym(ia)
      zmm(ia)=zm(ia)
      xm(ia)=x(ia)
      ym(ia)=y(ia)
1     zm(ia)=z(ia)
c
c     write Fourier coefficients for restart:
      if(mod(istep,iwrite).eq.0)then
       call rwfc(0,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux,
     1 wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux,
     3 tsuy,tcuy,tsuz,tcuz,
     4 wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0)
c
c      selected atoms
c      write the normalized power spectrum
       do 711 iw=1,nw
       w=wmin+dble(iw-1)*dw
       fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol
711    spec(iw)=fac*
     1(wsux(iw)*wsux(iw)+wsuy(iw)*wsuy(iw)+wsuz(iw)*wsuz(iw)
     2+wcux(iw)*wcux(iw)+wcuy(iw)*wcuy(iw)+wcuz(iw)*wcuz(iw))
c
c      convolute it with a smoothing function and write on disk
       call wrspec('WDraw.PRN',nw,wmin,dw,spec,tp,sol)
       call conv(res,nw,stem,spec,wmin,dw,tp,sol)
       call wrspec('WD.PRN',nw,wmin,dw,spec,tp,sol)

       do 8 iw=1,nw
       w=wmin+dble(iw-1)*dw
       fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol
8      spec(iw)=fac*
     1(wsux(iw)*wcmx(iw)+wsuy(iw)*wcmy(iw)+wsuz(iw)*wcmz(iw)
     2-wcux(iw)*wsmx(iw)-wcuy(iw)*wsmy(iw)-wcuz(iw)*wsmz(iw))
c
       call wrspec('WRraw.PRN',nw,wmin,dw,spec,tp,sol)
       call conv(res,nw,stem,spec,wmin,dw,tp,sol)
       call wrspec('WR.PRN',nw,wmin,dw,spec,tp,sol)

c      whole system response:
       do 9 iw=1,nw
       w=wmin+dble(iw-1)*dw
       fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol
9      spec(iw)=fac*
     1(tsux(iw)*tsux(iw)+tsuy(iw)*tsuy(iw)+tsuz(iw)*tsuz(iw)
     2+tcux(iw)*tcux(iw)+tcuy(iw)*tcuy(iw)+tcuz(iw)*tcuz(iw))
       call wrspec('WDTraw.PRN',nw,wmin,dw,spec,tp,sol)
       call conv(res,nw,stem,spec,wmin,dw,tp,sol)
       call wrspec('WDT.PRN',nw,wmin,dw,spec,tp,sol)
       do 11 iw=1,nw
       w=wmin+dble(iw-1)*dw
       fac=(w/(tp*sol)/(kelvin/0.69502d0))*w/tp**2/9.184d-3*dbb/tp/sol
11     spec(iw)=fac*
     1(tsux(iw)*tcmx(iw)+tsuy(iw)*tcmy(iw)+tsuz(iw)*tcmz(iw)
     2-tcux(iw)*tsmx(iw)-tcuy(iw)*tsmy(iw)-tcuz(iw)*tsmz(iw))
       call wrspec('WRTraw.PRN',nw,wmin,dw,spec,tp,sol)
       call conv(res,nw,stem,spec,wmin,dw,tp,sol)
       call wrspec('WRT.PRN',nw,wmin,dw,spec,tp,sol)
c
c      Raman-selected atoms only:
       if(ltti.or.lexp)then
        open(51,file='RAM.SCR',form='unformatted')
        do 1233 iw=1,nw
        w=wmin+dble(iw-1)*dw
c       temperature factor:*1d5(arbitrary factor)
        if(dabs(w).lt.1.0d-40)then
	 tf=0.0d0
	else
         tf=1.0d0/w/(1.0d0-exp(-w*1.44d0/298.0d0))*1.0d5
	endif
	if(lexp)then
	 tf=tf*w*1.0d-10
	else
	 tf=tf*1.0d6
	endif
        sal0=0.0d0
        sal1=0.0d0
        sag0=0.0d0
        sag1=0.0d0
        sa1=0.0d0
        do 121 jx=1,3
        sa1=sa1+wsat(1,jx,iw)*(wsaat(2,3,jx,iw)-wsaat(3,2,jx,iw))
     1         +wcat(1,jx,iw)*(wcaat(2,3,jx,iw)-wcaat(3,2,jx,iw))
     2         +wsat(2,jx,iw)*(wsaat(3,1,jx,iw)-wsaat(1,3,jx,iw))
     3         +wcat(2,jx,iw)*(wcaat(3,1,jx,iw)-wcaat(1,3,jx,iw))
     4         +wsat(3,jx,iw)*(wsaat(1,2,jx,iw)-wsaat(2,1,jx,iw))
     5         +wcat(3,jx,iw)*(wcaat(1,2,jx,iw)-wcaat(2,1,jx,iw))
        do 121 ix=1,3
        sal0=sal0+wsat(ix,ix,iw)*wsat(jx,jx,iw)
     1           +wcat(ix,ix,iw)*wcat(jx,jx,iw)
        sal1=sal1+wsat(ix,jx,iw)*wsat(ix,jx,iw)
     1           +wcat(ix,jx,iw)*wcat(ix,jx,iw)
c       coordinates
        sag0=sag0+wsat(ix,ix,iw)*wsgt(jx,jx,iw)
     1           +wcat(ix,ix,iw)*wcgt(jx,jx,iw)
        sag1=sag1+wsat(ix,jx,iw)*wsgt(ix,jx,iw)
     1           +wcat(ix,jx,iw)*wcgt(ix,jx,iw)
c       momenta:
c       sag0=sag0+wsat(ix,ix,iw)*wcgt(jx,jx,iw)
c    1           -wcat(ix,ix,iw)*wsgt(jx,jx,iw)
c       sag1=sag1+wsat(ix,jx,iw)*wcgt(ix,jx,iw)
c    1           -wcat(ix,jx,iw)*wsgt(ix,jx,iw)
121     continue
        sa1=sa1/3.0d0
c
c       IL+IR Backscattering
        s180=7.0d0*sal1+sal0
c       IL+IR Backscattering-xpol
        ydy=3.0d0*sal1-sal0
c       IL+IR Backscattering-ypol
        ydx=4.0d0*sal1+2.0d0*sal0
c       IR-IL Backscattering
        d180=8.0d0*(3.0d0*sag1-sag0+sa1)
1233    write(51)s180*tf,ydy*tf,ydx*tf,d180*tf
        rewind 51
        do 1234 iw=1,nw
        read(51)s180,ydy,ydx,d180
1234    spec(iw)=(w/sol/1000.0d0)**2/dble(nstep)**2*s180
        call wrspec('WRam180raw.PRN',nw,wmin,dw,spec,tp,sol)
        call conv(res,nw,stem,spec,wmin,dw,tp,sol)
        call wrspec('WRam180.PRN',nw,wmin,dw,spec,tp,sol)
        rewind 51
        do 1235 iw=1,nw
        read(51)s180,ydy,ydx,d180
1235    spec(iw)=(w/sol/1000.0d0)**2/dble(nstep)**2*ydy
        call wrspec('WRamyraw.PRN',nw,wmin,dw,spec,tp,sol)
        call conv(res,nw,stem,spec,wmin,dw,tp,sol)
        call wrspec('WRamy.PRN',nw,wmin,dw,spec,tp,sol)
        rewind 51
        do 1236 iw=1,nw
        read(51)s180,ydy,ydx,d180
1236    spec(iw)=(w/sol/1000.0d0)**2/dble(nstep)**2*ydx
        call wrspec('WRamxraw.PRN',nw,wmin,dw,spec,tp,sol)
        call conv(res,nw,stem,spec,wmin,dw,tp,sol)
        call wrspec('WRamx.PRN',nw,wmin,dw,spec,tp,sol)
        rewind 51
        do 1237 iw=1,nw
        read(51)s180,ydy,ydx,d180
1237    spec(iw)=(w/sol/1000.0d0)**2/dble(nstep)**2*d180
        call wrspec('WRoa180raw.PRN',nw,wmin,dw,spec,tp,sol)
        call conv(res,nw,stem,spec,wmin,dw,tp,sol)
        call wrspec('WRoa180.PRN',nw,wmin,dw,spec,tp,sol)
        close(51)
c
       endif
      endif

      endif
c     ########################################################
      end do
      if(lfilterdiag)close(66)
      if(lexsp)close(67)
c
c     perform any final tasks before program exit
c
      call final
      end



      subroutine conv(res,nw,stem,spec,wmin,dw,tp,sol)
      implicit none
      real*8 res,stem(*),spec(*),wmin,dw,fac,tp,sol,w,wp
      integer*4 nw,iw,iwp
      if(res.gt.0)then
       do 72 iw=1,nw
72     stem(iw)=spec(iw)
       do 73 iw=1,nw
       w=(wmin+dble(iw-1)*dw)/tp/sol
       spec(iw)=0.0d0
       do 73 iwp=1,nw
       wp=(wmin+dble(iwp-1)*dw)/tp/sol
       fac=((wp-w)/res)**2
73     if(fac.lt.20.0d0)spec(iw)=spec(iw)+stem(iwp)*exp(-fac)
      endif
      return
      end

      subroutine wrspec(sn,nw,wmin,dw,spec,tp,sol)
      implicit none
      character*(*) sn
      real*8 wmin,dw,spec(*),w,tp,sol
      integer*4 nw,iw
      open(8,file=sn)
      do 1 iw=1,nw
      w=wmin+dble(iw-1)*dw
1     write(8,805)w/tp/sol,spec(iw)
805   format(f12.2,f25.12)
      close(8)
      return
      end

      subroutine update(wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,
     1wsux,wcux,wsuy,wcuy,wsuz,wcuz,mx,my,mz,ux,uy,uz,t,sc)
      implicit none
      integer*4 nw,iw
      real*8 wmin,w,mx,my,mz,ux,uy,uz,t,ss,cc,sc,dw,
     2wsux(*),wsuy(*),wsuz(*),wsmx(*),wsmy(*),wsmz(*),
     3wcux(*),wcuy(*),wcuz(*),wcmx(*),wcmy(*),wcmz(*)
      w=wmin-dw
      do 1 iw=1,nw
      w=w+dw
      ss=sin(w*t)
      cc=cos(w*t)
      wsmx(iw)=wsmx(iw)*sc+mx*ss
      wcmx(iw)=wcmx(iw)*sc+mx*cc
      wsmy(iw)=wsmy(iw)*sc+my*ss
      wcmy(iw)=wcmy(iw)*sc+my*cc
      wsmz(iw)=wsmz(iw)*sc+mz*ss
      wcmz(iw)=wcmz(iw)*sc+mz*cc
      wsux(iw)=wsux(iw)*sc+ux*ss
      wcux(iw)=wcux(iw)*sc+ux*cc
      wsuy(iw)=wsuy(iw)*sc+uy*ss
      wcuy(iw)=wcuy(iw)*sc+uy*cc
      wsuz(iw)=wsuz(iw)*sc+uz*ss
1     wcuz(iw)=wcuz(iw)*sc+uz*cc
      return
      end

      subroutine update2(wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,
     1wsux,wcux,wsuy,wcuy,wsuz,wcuz,mx,my,mz,ux,uy,uz,t,sc)
      implicit none
      integer*4 nw,iw
      real*8 wmin,w,mx,my,mz,ux,uy,uz,t,ss,cc,sc,dw,tp,sol,
     2wsux(*),wsuy(*),wsuz(*),wsmx(*),wsmy(*),wsmz(*),
     3wcux(*),wcuy(*),wcuz(*),wcmx(*),wcmy(*),wcmz(*),wau
      w=wmin-dw
      tp=2.0d0*4.0d0*datan(1.0d0)
      sol=3.0d10
      do 1 iw=1,nw
      w=w+dw
      wau=w/tp/sol/219470.0d0
      if(abs(w).gt.1.0d-10)then
       ss=-cos(w*t)/wau
       cc= sin(w*t)/wau
       wsmx(iw)=wsmx(iw)*sc+mx*ss
       wcmx(iw)=wcmx(iw)*sc+mx*cc
       wsmy(iw)=wsmy(iw)*sc+my*ss
       wcmy(iw)=wcmy(iw)*sc+my*cc
       wsmz(iw)=wsmz(iw)*sc+mz*ss
       wcmz(iw)=wcmz(iw)*sc+mz*cc
       wsux(iw)=wsux(iw)*sc+ux*ss
       wcux(iw)=wcux(iw)*sc+ux*cc
       wsuy(iw)=wsuy(iw)*sc+uy*ss
       wcuy(iw)=wcuy(iw)*sc+uy*cc
       wsuz(iw)=wsuz(iw)*sc+uz*ss
       wcuz(iw)=wcuz(iw)*sc+uz*cc
      endif
1     continue
      return
      end

      subroutine update3(wmin,dw,nw,wsat,wcat,wsgt,wcgt,
     1wsaat,wcaat,at,gt,aat,t,sc,MX3)
      implicit none
      integer*4 nw,iw,ix,jx,kx,MX3
      real*8 wmin,w,at(3,3),gt(3,3),aat(3,3,3),t,ss,cc,sc,dw,tp,sol,
     2wsat(3,3,MX3),wcat(3,3,MX3),wsgt(3,3,MX3),wcgt(3,3,MX3),
     2wsaat(3,3,3,MX3),wcaat(3,3,3,MX3),wau
      tp=2.0d0*4.0d0*datan(1.0d0)
      sol=3.0d10
      w=wmin-dw
      do 1 iw=1,nw
      w=w+dw
      wau=w/tp/sol/219470.0d0
      if(abs(w).gt.1.0d-10)then
       ss=-cos(w*t)/wau
       cc= sin(w*t)/wau
       do 2 ix=1,3
       do 2 jx=1,3
       wsat(ix,jx,iw)=wsat(ix,jx,iw)*sc+at(ix,jx)*ss
       wcat(ix,jx,iw)=wcat(ix,jx,iw)*sc+at(ix,jx)*cc
       wsgt(ix,jx,iw)=wsgt(ix,jx,iw)*sc+gt(ix,jx)*ss
       wcgt(ix,jx,iw)=wcgt(ix,jx,iw)*sc+gt(ix,jx)*cc
       do 2 kx=1,3
       wsaat(ix,jx,kx,iw)=wsaat(ix,jx,kx,iw)*sc+aat(ix,jx,kx)*ss
2      wcaat(ix,jx,kx,iw)=wcaat(ix,jx,kx,iw)*sc+aat(ix,jx,kx)*cc
      endif
1     continue
      return
      end

      subroutine ifc(nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux,
     1wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux,
     3tsuy,tcuy,tsuz,tcuz,
     4 wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0)
      implicit none
      integer iw,nw,ix,jx,kx,nw0
      real*8 wsmx(*),wcmx(*),wsmy(*),wcmy(*),wsmz(*),wcmz(*),
     1wsux(*),wcux(*),wsuy(*),wcuy(*),wsuz(*),wcuz(*),tsmx(*),
     2tcmx(*),tsmy(*),tcmy(*),tsmz(*),tcmz(*),tsux(*),tcux(*),
     3tsuy(*),tcuy(*),tsuz(*),tcuz(*),
     2wsat(3,3,nw0),wcat(3,3,nw0),wsgt(3,3,nw0),wcgt(3,3,nw0),
     2wsaat(3,3,3,nw0),wcaat(3,3,3,nw0)
      do 1 iw=1,nw
      wsmx(iw)=0.0d0
      wcmx(iw)=0.0d0
      wsmy(iw)=0.0d0
      wcmy(iw)=0.0d0
      wsmz(iw)=0.0d0
      wcmz(iw)=0.0d0
      wsux(iw)=0.0d0
      wcux(iw)=0.0d0
      wsuy(iw)=0.0d0
      wcuy(iw)=0.0d0
      wsuz(iw)=0.0d0
      wcuz(iw)=0.0d0
      tsmx(iw)=0.0d0
      tcmx(iw)=0.0d0
      tsmy(iw)=0.0d0
      tcmy(iw)=0.0d0
      tsmz(iw)=0.0d0
      tcmz(iw)=0.0d0
      tsux(iw)=0.0d0
      tcux(iw)=0.0d0
      tsuy(iw)=0.0d0
      tcuy(iw)=0.0d0
      tsuz(iw)=0.0d0
      tcuz(iw)=0.0d0
      do 1 ix=1,3
      do 1 jx=1,3
      wsat(ix,jx,iw)=0.0d0
      wcat(ix,jx,iw)=0.0d0
      wsgt(ix,jx,iw)=0.0d0
      wcgt(ix,jx,iw)=0.0d0
      do 1 kx=1,3
      wcaat(ix,jx,kx,iw)=0.0d0
1     wsaat(ix,jx,kx,iw)=0.0d0
      return
      end

      subroutine rwfc(ic,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,wsux,wcux,
     1wsuy,wcuy,wsuz,wcuz,tsmx,tcmx,tsmy,tcmy,tsmz,tcmz,tsux,tcux,
     3tsuy,tcuy,tsuz,tcuz,
     4 wsat,wcat,wsgt,wcgt,wsaat,wcaat,nw0)
      implicit none
      integer iw,nw,ic,nw0,ix,jx,kx
      real*8 wsmx(*),wcmx(*),wsmy(*),wcmy(*),wsmz(*),wcmz(*),
     1wsux(*),wcux(*),wsuy(*),wcuy(*),wsuz(*),wcuz(*),tsmx(*),
     2tcmx(*),tsmy(*),tcmy(*),tsmz(*),tcmz(*),tsux(*),tcux(*),
     3tsuy(*),tcuy(*),tsuz(*),tcuz(*),
     2wsat(3,3,nw0),wcat(3,3,nw0),wsgt(3,3,nw0),wcgt(3,3,nw0),
     2wsaat(3,3,3,nw0),wcaat(3,3,3,nw0)
      open(9,file='FIT.SCR',form='unformatted')
      if(ic.eq.0)then
      write(9)
     1(wsmx(iw),iw=1,nw),(wsmy(iw),iw=1,nw),(wsmz(iw),iw=1,nw),
     1(wcmx(iw),iw=1,nw),(wcmy(iw),iw=1,nw),(wcmz(iw),iw=1,nw),
     2(wsux(iw),iw=1,nw),(wsuy(iw),iw=1,nw),(wsuz(iw),iw=1,nw),
     2(wcux(iw),iw=1,nw),(wcuy(iw),iw=1,nw),(wcuz(iw),iw=1,nw),
     3(tsmx(iw),iw=1,nw),(tsmy(iw),iw=1,nw),(tsmz(iw),iw=1,nw),
     3(tcmx(iw),iw=1,nw),(tcmy(iw),iw=1,nw),(tcmz(iw),iw=1,nw),
     4(tsux(iw),iw=1,nw),(tsuy(iw),iw=1,nw),(tsuz(iw),iw=1,nw),
     4(tcux(iw),iw=1,nw),(tcuy(iw),iw=1,nw),(tcuz(iw),iw=1,nw),
     4(((wsat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4(((wcat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4(((wsgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4(((wcgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4((((wsaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw),
     4((((wcaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw)
      else
      read(9)
     1(wsmx(iw),iw=1,nw),(wsmy(iw),iw=1,nw),(wsmz(iw),iw=1,nw),
     1(wcmx(iw),iw=1,nw),(wcmy(iw),iw=1,nw),(wcmz(iw),iw=1,nw),
     2(wsux(iw),iw=1,nw),(wsuy(iw),iw=1,nw),(wsuz(iw),iw=1,nw),
     2(wcux(iw),iw=1,nw),(wcuy(iw),iw=1,nw),(wcuz(iw),iw=1,nw),
     3(tsmx(iw),iw=1,nw),(tsmy(iw),iw=1,nw),(tsmz(iw),iw=1,nw),
     3(tcmx(iw),iw=1,nw),(tcmy(iw),iw=1,nw),(tcmz(iw),iw=1,nw),
     4(tsux(iw),iw=1,nw),(tsuy(iw),iw=1,nw),(tsuz(iw),iw=1,nw),
     4(tcux(iw),iw=1,nw),(tcuy(iw),iw=1,nw),(tcuz(iw),iw=1,nw),
     4(((wsat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4(((wcat(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4(((wsgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4(((wcgt(ix,jx,iw),ix=1,3),jx=1,3),iw=1,nw),
     4((((wsaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw),
     4((((wcaat(ix,jx,kx,iw),ix=1,3),jx=1,3),kx=1,3),iw=1,nw)
      endif
      close(9)
      return
      end

      SUBROUTINE READFF(MX3,N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(MX3,N)
      open(20,file='FILE.FC')
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
C
c     CONST=4.359828d0/0.5291772d0/0.5291772d0
      CONST=627.5d0/0.5291772d0/0.5291772d0
      DO 3 I=1,N
      DO 3 J=I,N
3     FCAR(J,I)=FCAR(J,I)*CONST
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      close(20)
      WRITE(6,*)' Cartesian FF read in from FILE.FC'
      RETURN
      END


      SUBROUTINE projectff(MX3,NA,FCAR,xeq,wfmax,wsup,mass,e,
     1fnew,smat,Q,ZM,ZIM,TEM)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(MX3,*),xeq(*),ZIM(*),ZM(*),
     1FNEW(MX3,MX3),SMAT(MX3,MX3),Q(*),TEM(MX3,MX3)
      real*8 mass(*),e(*)
C
      NOAT=NA/3

      BM=0.0d0
      DO 102 IA=1,NOAT
      ZIM(IA)=mass(IA)
      DO 220 J=1,3
220   ZM(3*IA-3+J)=1.0d0/SQRT(ZIM(IA))
102   BM=BM+ZIM(IA)
      WRITE(6,778)BM
778   FORMAT(' TOTAL MASS  ',F20.6)
C
      CALL PROJECT(FCAR,MX3,NA,FNEW,TEM,xeq)
      write(6,*)' Translation-rotations projected'
C
      CALL FAST(NA,FCAR,FNEW,ZM,MX3)
      IMOD=2
      CALL TRED12(MX3,NA,FNEW,SMAT,Q,IMOD,IERR,e)
      IF(IERR.NE.0)STOP
      WRITE(6,*)' Diagonalization OK'

      DO 310 I=1,NA
      T=Q(I)
      SIGN=1.0d0
      IF(T.LT.0.0d0)SIGN=-1.0d0
      T=SIGN*SQRT(ABS(T))*1302.828d0
  310 Q(I)=T
C
c     set low modes to required constant:
      WRITE(*,5555)
5555  FORMAT(1X,'************  NEW FREQUENCIES ******************')
      DO 311 J=NA,1,-1
      if(Q(NA-J+1).lt.wfmax)then
       WRITE(*,6466)J,NA-J+1,Q(NA-J+1),wsup
       Q(NA-J+1)=(wsup/1302.828d0)**2
6466   FORMAT(I6,' (',I4,') ...... ',F8.2,' changed to ',f8.2)
      else
       WRITE(*,6666)J,NA-J+1,Q(NA-J+1)
6666   FORMAT(I6,' (',I4,') ...... ',F8.2,' (unchanged)')
       Q(NA-J+1)=(Q(NA-J+1)/1302.828d0)**2
      endif
311   continue
      WRITE(*,5551)
5551  FORMAT(1X,'************************************************')
C
c     project back:
      do 1 i=1,NA
      do 1 j=1,na
      FNEW(i,j)=0.0d0
      do 1 l=1,na
1     FNEW(i,j)=FNEW(i,j)+SMAT(i,l)*Q(l)*SMAT(j,l)

      CALL FASTm1(NA,FCAR,FNEW,ZM,MX3)

      write(6,*)' New force field made'
      return
      END
C
C
      SUBROUTINE FAST(NA,FCAR,FNEW,ZM,MX3)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      DIMENSION FCAR(MX3,NA),FNEW(MX3,NA),ZM(NA)
      CONST=627.5d0/4.359828d0
      DO 214 I=1,NA
      DO 214 J=I,NA
      FNEW(I,J)=FCAR(I,J)*ZM(J)*ZM(I)/CONST
214   FNEW(J,I)=FNEW(I,J)
      RETURN
      END
C
      SUBROUTINE FASTm1(NA,FCAR,FNEW,ZM,MX3)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      DIMENSION FCAR(MX3,NA),FNEW(MX3,NA),ZM(NA)
      CONST=627.5d0/4.359828d0
      DO 214 I=1,NA
      DO 214 J=I,NA
      FCAR(I,J)=FNEW(I,J)/ZM(J)/ZM(I)*CONST
214   FCAR(J,I)=FCAR(I,J)
      RETURN
      END
C
      SUBROUTINE PROJECT(F,N0,N,A,TEM,xeq)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION F(N0,N),A(N0,N0),TEM(N0,N0),xeq(*)
      CALL DOMA(A,N0,N,xeq,TEM)
      DO 1 I=1,N
      DO 1 J=1,N
      TEM(I,J)=0.0d0
      DO 1 II=1,N
1     TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J)
      DO 2 I=1,N
      DO 2 J=1,N
      F(I,J)=0.0d0
      DO 2 II=1,N
2     F(I,J)=F(I,J)+TEM(I,II)*A(J,II)
      RETURN
      END

      SUBROUTINE DOMA(A,NDIM,N,xeq,TEM)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(NDIM,N),TEM(NDIM,*),xeq(*)
C
      AMACH=1.0d-14
      NAT=N/3
      DO 12 I=1,N
      DO 12 J=1,6
12    A(I,J)=0.0d0
      DO 1 I=1,NAT
      A(3*I-2,1)= 1.0d0
      A(3*I-1,2)= 1.00d0
      A(3*I  ,3)= 1.00d0
      A(3*I-2,4)=-xeq(3*(I-1)+2)
      A(3*I-1,4)= xeq(3*(I-1)+1)
      A(3*I-1,5)=-xeq(3*(I-1)+3)
      A(3*I  ,5)= xeq(3*(I-1)+2)
      A(3*I-2,6)= xeq(3*(I-1)+3)
1     A(3*I  ,6)=-xeq(3*(I-1)+1)
      N6=6
      CALL ORT(A,NDIM,N6,NAT)
      ICOL=1
5     IF(ICOL.GT.N6)GOTO 8
      S=0.0d0
      DO 6 J=1,N
6     S=S+A(J,ICOL)**2
      IF(S.LT.1000.0d0*AMACH)THEN
       DO 7 J=1,N
7      A(J,ICOL)=A(J,N6)
       N6=N6-1
       GOTO 5
      ENDIF
      ICOL=ICOL+1
      GOTO 5
8     CONTINUE
      WRITE(6,*)N6,' coordinates projected'
      DO 2 I=1,N
      DO 2 J=1,N
      TEM(I,J)=0.0d0
      DO 2 II=1,N6
2     TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II)
      DO 3 I=1,N
      DO 4 J=1,N
4     A(I,J)=-TEM(I,J)
3     A(I,I)=A(I,I)+1.0d0
      RETURN
      END
C     ============================================================ 
      SUBROUTINE ORT(A,NDIM,N6,NAT)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(NDIM,N6)
      AMACH=1.0d-14
      NAT3=3*NAT
      C=0.0d0
      DO 1 I=1,NAT3
1     C=C+A(I,1)**2
      IF(C.GT.AMACH)C=1.0d0/SQRT(C)
      DO 2 I=1,NAT3
2     A(I,1)=A(I,1)*C
      DO 3 IC=2,N6
C     Column IC, orthonormalize to previous columns:
      DO 3 ICP=1,IC-1
      SP=0.0d0
      DO 5 J=1,NAT3
5     SP=SP+A(J,IC)*A(J,ICP)
C     Subtract projection to ICP:
      DO 6 J=1,NAT3
6     A(J,IC)=A(J,IC)-SP*A(J,ICP)
C     Normalize :
      SP=0.0d0
      DO 7 J=1,NAT3
7     SP=SP+A(J,IC)**2
      IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP)
      DO 3 J=1,NAT3
3     A(J,IC)=A(J,IC)*SP
      RETURN
      END
C
      SUBROUTINE READTEN(nat0ai,P,A)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(nat0ai,3,3),A(nat0ai,3,3)
      open(15,file='FILE.TEN')
      READ (15,*) NAT
      DO 10 L=1,NAT
      DO 10 J=1,3                                                     
10    READ (15,*) (P(L,J,I),I=1,3)
      DO 220 L=1,NAT                                                 
      DO 220 J=1,3                                                    
220   READ (15,*) (A(L,J,I),I=1,3)
      close(15)
      write(6,*)'FILE.TEN read'
      RETURN
      END

      SUBROUTINE READTTT(MX,ALPHA,A,G,NAT)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3)
      READ(15,*)
      READ(15,*)NAT
      READ(15,*)
      READ(15,*)
      DO 1 I=1,3
      READ(15,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IIND=3*(L-1)+IX
1     READ(15,*)LDUM,IXDUM,(ALPHA(IIND,I,J),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 2 I=1,3
      READ(15,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
2     READ(15,*)LDUM,IXDUM,(G(IIND,I,J),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(15,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IIND=3*(L-1)+IX
3     READ(15,*)LDUM,IXDUM,(A(IIND,I,J,K),K=1,3)
      RETURN
      END

      SUBROUTINE READTENI(N0,P,A,NAT)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (nat0=1000)
      DIMENSION P(N0,3,3),A(N0,3,3),AI(nat0,3,3),AJ(nat0,3,3)
      READ (15,*) NOAT
      if(NOAT.ne.NAT)call report('NAT<>NOAT')
      DO 10 L=1,NAT
      DO 10 J=1,3                                                     
10    READ (15,*) (P(L,J,I),I=1,3)
      DO 220 L=1,NOAT                                                 
      DO 220 J=1,3                                                    
220   READ (15,*) (AI(L,J,I),I=1,3)
      DO 230 L=1,NOAT                                                 
      DO 230 J=1,3                                                    
230   READ (15,*) (AJ(L,J,I),I=1,3)                                   
      DO 100 L=1,NOAT                                                 
      DO 100 J=1,3                                                    
100   READ (15,*)
C     
C     Total axial tensor = electronic + nuclear:
      DO 4 L=1,NAT
      DO 4 I=1,3
      DO 4 J=1,3
4     A(L,I,J)=AI(L,I,J)+AJ(L,I,J)
C
      RETURN
      END

      subroutine readatms(io,NAT,iii,MX)
      integer*4 MX,NAT,iii(4,MX),ia,io,idumm
      do 1 ia=1,NAT
1     read(io,*)idumm,(iii(ii,ia),ii=1,4)
      return
      end

      subroutine dtm(ii,NAT,u,MX,x,y,z)
      implicit none
      integer*4 MX,NAT,ii(4,MX),ia,i1,i2,j1,j2,ix,jx
      real*8 u(3,3,MX),x(MX),v1(3),v2(3),v3(3),y(MX),z(MX)
c
      do 1 ia=1,NAT
      i1=ii(1,ia)
      i2=ii(2,ia)
      j1=ii(3,ia)
      j2=ii(4,ia)
      if(i1.gt.0.and.i2.gt.0.and.j1.gt.0.and.j2.gt.0)then
       v1(1)=x(i2)-x(i1)
       v2(1)=x(j2)-x(j1)
       v1(2)=y(i2)-y(i1)
       v2(2)=y(j2)-y(j1)
       v1(3)=z(i2)-z(i1)
       v2(3)=z(j2)-z(j1)
       v3(1)=v1(2)*v2(3)-v1(3)*v2(2)
       v3(2)=v1(3)*v2(1)-v1(1)*v2(3)
       v3(3)=v1(1)*v2(2)-v1(2)*v2(1)
       call norm(v1)
       call norm(v3)
       v2(1)=v3(2)*v1(3)-v3(3)*v1(2)
       v2(2)=v3(3)*v1(1)-v3(1)*v1(3)
       v2(3)=v3(1)*v1(2)-v3(2)*v1(1)
       do 13 ix=1,3
       u(1,ix,ia)=v1(ix)
       u(2,ix,ia)=v2(ix)
13     u(3,ix,ia)=v3(ix)
      else
       do 14 ix=1,3
       do 14 jx=1,3
14     u(jx,ix,ia)=0.0d0
      endif
1     continue
      return
      end

      SUBROUTINE TTTT(MX3,ALPHA,A,G,NAT,ICO,xx,y,z,X,MX)
C     Transforms A and G to local origins (ICO=1) 
c     or common origin (ICO=2);
C     sets A and G to zero for ICO=3.
C     R supposed to be passed in in A
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(MX3,3,3),G(MX3,3,3),A(MX3,3,3,3),xx(*),X(3,MX),
     1y(*),z(*)
C
C     transfer into atomic coordinates:
      DO 5 J=1,NAT
      X(1,J)=xx(J)/0.529177D0
      X(2,J)= y(J)/0.529177D0
5     X(3,J)= z(J)/0.529177D0
C
      IF(ICO.EQ.1)SIGN=1.0d0
      IF(ICO.EQ.2)SIGN=-1.0d0
      IF(ICO.EQ.3)SIGN=0.0d0
      W=1.0d0
C     supposedly the static limit G/w is calculated
C
      DO 2 IA=1,3
      DO 2 IB=1,3
      DO 2 L=1,NAT
      DO 2 IE=1,3
      IIND=3*(L-1)+IE
      SUM=0.0d0
      DO 1 IG=1,3
      DO 1 ID=1,3
1     SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID)
      G(IIND,IA,IB)=G(IIND,IA,IB)+SUM
      IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0
2     CONTINUE
c     IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins'
c     IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin'
c     IF(ICO.EQ.3)WRITE(6,*)' G set to zero'
C
      DO 3 IA=1,3
      DO 3 IB=1,3
      DO 3 IC=1,3
      DO 3 L=1,NAT
      DO 3 IE=1,3
      IIND=3*(L-1)+IE
      SUM=0.0d0
      IF(IB.EQ.IC)THEN
       DO 4 ID=1,3
4      SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN
      ENDIF
      A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM
     1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB))
      IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0
3     CONTINUE
c     IF(ICO.EQ.1)WRITE(6,*)' A transformed into atomic origins'
c     IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin'
c     IF(ICO.EQ.3)WRITE(6,*)' A set to zero'
      RETURN
      END

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      FUNCTION EPS(I,J,K)
      REAL*8 EPS
      INTEGER*4 I,J,K
      EPS=0.0d0
      IF(I.EQ.1)then
       if(J.EQ.2.AND.K.EQ.3)then
        EPS= 1.0d0
       else
        IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0d0
       endif
      else
       IF(I.EQ.2)then
        IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0d0
        IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0d0
       else
        IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0d0
        IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0d0
       endif
      endif
      RETURN                                                          
      END                                                             

      subroutine norm(v)
      implicit none
      real*8 v(*),n
      n=dsqrt(v(1)**2+v(2)**2+v(3)**2)
      v(1)=v(1)/n
      v(2)=v(2)/n
      v(3)=v(3)/n
      return
      end

      subroutine updt(na1,na2,wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,
     1wcmz,wsux,wcux,wsuy,wcuy,wsuz,wcuz,t,sc,xm,ym,zm,xmm,ymm,zmm,
     2bohr,dtau,x,y,z,A,P,nat0ai,ltei,lten,pchgft,ux)
      implicit none
      logical ltei,lten
      integer*4 na1,na2,ia,nw,nat0ai
      real*8 wmin,dw,mx,my,mz,ux,uy,uz,t,sc,bohr,dtau,x(*),y(*),z(*),
     2xm(*),ym(*),zm(*),xmm(*),ymm(*),zmm(*),cf1,cf2,vx,vy,vz,ax,ay,az,
     3P(nat0ai,3,3),A(nat0ai,3,3),qia,pchgft(*),fac,
     2wsux(*),wsuy(*),wsuz(*),wsmx(*),wsmy(*),wsmz(*),
     3wcux(*),wcuy(*),wcuz(*),wcmx(*),wcmy(*),wcmz(*)
c
      mx=0.0d0
      my=0.0d0
      mz=0.0d0
      ux=0.0d0
      uy=0.0d0
      uz=0.0d0
      if(lten.or.ltei)then
c      constants changing vx,ax into atomic units
       cf1=1.0d0/bohr/dtau
       cf2=1.0d0/bohr/dtau/dtau
       do 201 ia=na1,na2
       vx=(x(ia)-xm(ia))*cf1
       vy=(y(ia)-ym(ia))*cf1
       vz=(z(ia)-zm(ia))*cf1
       ax=(x(ia)+xmm(ia)-xm(ia)*2.0d0)*cf2
       ay=(y(ia)+ymm(ia)-ym(ia)*2.0d0)*cf2
       az=(z(ia)+zmm(ia)-zm(ia)*2.0d0)*cf2
       mx=mx+ax*A(ia,1,1)+ay*A(ia,2,1)+az*A(ia,3,1)
       my=my+ax*A(ia,1,2)+ay*A(ia,2,2)+az*A(ia,3,2)
       mz=mz+ax*A(ia,1,3)+ay*A(ia,2,3)+az*A(ia,3,3)
       ux=ux+vx*P(ia,1,1)+vy*P(ia,2,1)+vz*P(ia,3,1)
       uy=uy+vx*P(ia,1,2)+vy*P(ia,2,2)+vz*P(ia,3,2)
201    uz=uz+vx*P(ia,1,3)+vy*P(ia,2,3)+vz*P(ia,3,3)
       call update2(wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,
     1 wsux,wcux,wsuy,wcuy,wsuz,wcuz,mx,my,mz,ux,uy,uz,t,sc)
      else
       cf1=1.0d0/2.0d0/dtau/bohr/bohr
       do 2 ia=na1,na2
       qia=pchgft(ia)
       fac=qia*cf1
       mx=mx+(y(ia)*(z(ia)-zm(ia))-z(ia)*(y(ia)-ym(ia)))*fac
       my=my+(z(ia)*(x(ia)-xm(ia))-x(ia)*(z(ia)-zm(ia)))*fac
       mz=mz+(x(ia)*(y(ia)-ym(ia))-y(ia)*(x(ia)-xm(ia)))*fac
       fac=qia/bohr
       ux=ux+x(ia)*fac
       uy=uy+y(ia)*fac
2      uz=uz+z(ia)*fac
       call update(wmin,dw,nw,wsmx,wcmx,wsmy,wcmy,wsmz,wcmz,
     1 wsux,wcux,wsuy,wcuy,wsuz,wcuz,mx,my,mz,ux,uy,uz,t,sc)
      endif
      return
      end


      include 'cctn_tinker1.f'