**********************************************************************
*                                                                    *
*     REAXFF Reactive force field program                            *
*                                                                    *
*     Developed and written by Adri van Duin, duin@wag.caltech.edu   *
*                                                                    *
*     Copyright (c) 2001-2010 California Institute of Technology     *
*                                                                    *
*     This is an open-source program. Feel free to modify its        *
*     contents. Please keep me informed of any useful modification   *
*     or addition that you made. Please do not distribute this       *
*     program to others; if people are interested in obtaining       *
*     a copy of this program let them contact me first.              *
*                                                                    *
**********************************************************************
**********************************************************************

      subroutine ffinpt

**********************************************************************
#include "cbka.blk"
#include "cbkboncor.blk"
#include "cbkconst.blk"
#include "cbkcovbon.blk"
#include "cbkff.blk"
#include "cbkfftorang.blk"
#include "cbknonbon.blk"
#include "cbksrthb.blk"
#include "cbktorsion.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "valang.blk"
#include "cbksrtbon1.blk"
#include "cbkchb.blk"
      dimension rcore2(nsort),ecore2(nsort),acore2(nsort)
**********************************************************************
*                                                                    *
*     Read in force field                                            *
*                                                                    *
**********************************************************************
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In ffinpt'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      open (4,file='ffield.reax',status='old')
      rewind (4)
      iline=0
      read (4,'(a40)',end=990,err=990)qffield
      iline=iline+1
**********************************************************************
*                                                                    *
*     Read in general force field parameters                         *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990)npar
      iline=iline+1
      do i1=1,npar
      read (4,1300,end=990,err=990)vpar(i1)
      iline=iline+1
      end do
      cutoff=0.01*vpar(30)
      swa=vpar(12)
      if (abs(swa).gt.0.01) write (*,*)
     $'Warning: non-zero value for lower Taper-radius cutoff'
      swb=vpar(13)
      if (swb.lt.zero) stop 
     $'Negative value for upper Taper-radius cutoff'
      if (swb.lt.5.0) write (*,*)
     $'Warning: very low value for upper Taper-radius cutoff:',swb
**********************************************************************
*                                                                    *
*     Read in atom type data                                         *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990) nso
      iline=iline+1
      read (4,*,end=990,err=990)
      iline=iline+1
      read (4,*,end=990,err=990)
      iline=iline+1
      read (4,*,end=990,err=990)
      iline=iline+1
      if (nso.gt.nsort) stop 'Maximum number of atom types exceeded'
      do i1=1,nso
      read (4,1200,end=990,err=990)qas(i1),rat(i1),aval(i1),amas(i1),
     $rvdw(i1),eps(i1),gam(i1),rapt(i1),stlp(i1)
      iline=iline+1
      read (4,1250,end=990,err=990)alf(i1),vop(i1),valf(i1),
     $valp1(i1),valp2(i1),chi(i1),eta(i1),vnphb
      iline=iline+1
      read (4,1250,end=990,err=990)vnq(i1),vlp1(i1),vincr(i1),
     $bo131(i1),bo132(i1),bo133(i1),sigqeq(i1),default
      iline=iline+1
      read (4,1250,end=990,err=990)vovun(i1),vval1(i1),vrom,
     $vval3(i1),vval4(i1),rcore2(i1),ecore2(i1),acore2(i1)
      iline=iline+1
      idef(i1)=int(default)
      nphb(i1)=int(vnphb)
      end do
**********************************************************************
*                                                                    *
*     Calculate van der Waals and Coulomb pair-parameters            *
*                                                                    *
**********************************************************************
      do i1=1,nso
      do i2=1,nso
      rcore(i1,i2)=sqrt(rcore2(i1)*rcore2(i2))
      ecore(i1,i2)=sqrt(ecore2(i1)*ecore2(i2))
      acore(i1,i2)=sqrt(acore2(i1)*acore2(i2))
      p1co(i1,i2)=sqrt(4.0*rvdw(i1)*rvdw(i2))
      p2co(i1,i2)=sqrt(eps(i1)*eps(i2))
      p3co(i1,i2)=sqrt(alf(i1)*alf(i2))
      gamwh=sqrt(vop(i1)*vop(i2))
      gamwco(i1,i2)=1.0/gamwh**vpar(29)
      gamch=sqrt(gam(i1)*gam(i2))
      gamcco(i1,i2)=1.0/gamch**3
      rob1(i1,i2)=0.50*(rat(i1)+rat(i2))
      rob2(i1,i2)=0.50*(rapt(i1)+rapt(i2))
      rob3(i1,i2)=0.50*(vnq(i1)+vnq(i2))
      end do
      end do
**********************************************************************
*                                                                    *
*     Read in bond type data                                         *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990)nboty
      iline=iline+1
      read (4,*,end=990,err=990)
      iline=iline+1
      if (2*nboty.gt.nbotym) stop 'Maximum nr. of bond types exceeded'
      ih=0
      do i1=1,nboty
      ih=ih+1
      read (4,1400,end=990,err=990)nbs(ih,1),nbs(ih,2),de1(ih),
     $de2(ih),de3(ih),psi(ih),pdo(ih),v13cor(ih),popi(ih),vover(ih)
      iline=iline+1
      read (4,1450,end=990,err=990)psp(ih),pdp(ih),ptp(ih),
     $bom(ih),bop1(ih),bop2(ih),ovc(ih),vuncor(ih)
      iline=iline+1
      if (nbs(ih,1).ne.nbs(ih,2)) then
      ih=ih+1
      nbs(ih,1)=nbs(ih-1,2)
      nbs(ih,2)=nbs(ih-1,1)
      de1(ih)=de1(ih-1)
      de2(ih)=de2(ih-1)
      de3(ih)=de3(ih-1)
      psi(ih)=psi(ih-1)
      pdo(ih)=pdo(ih-1)
      v13cor(ih)=v13cor(ih-1)
      vover(ih)=vover(ih-1)
      psp(ih)=psp(ih-1)
      pdp(ih)=pdp(ih-1)
      ptp(ih)=ptp(ih-1)
      bop1(ih)=bop1(ih-1)
      bop2(ih)=bop2(ih-1)
*     bop3(ih)=bop3(ih-1)
*     bop4(ih)=bop4(ih-1)
      bom(ih)=bom(ih-1)
      popi(ih)=popi(ih-1)
      ovc(ih)=ovc(ih-1)
      end if
      end do
      nboty2=ih
**********************************************************************
*                                                                    *
*     Read in off-diagonal parameters                                *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990)nodmty
      iline=iline+1
      if (nodmty.gt.nodmtym)
     $stop 'Maximum nr. of off-diagonal Morse types exceeded'
      ih=0
      do i1=1,nodmty
      ih=ih+1
      read (4,1400,end=990,err=990)nodm1,nodm2,deodmh,rodmh,godmh,
     $rsig,rpi,rpi2
      iline=iline+1
      if (rsig.gt.zero) rob1(nodm1,nodm2)=rsig
      if (rsig.gt.zero) rob1(nodm2,nodm1)=rsig
      if (rpi.gt.zero) rob2(nodm1,nodm2)=rpi
      if (rpi.gt.zero) rob2(nodm2,nodm1)=rpi
      if (rpi2.gt.zero) rob3(nodm1,nodm2)=rpi2
      if (rpi2.gt.zero) rob3(nodm2,nodm1)=rpi2
      if (rodmh.gt.zero) p1co(nodm1,nodm2)=2.0*rodmh
      if (rodmh.gt.zero) p1co(nodm2,nodm1)=2.0*rodmh
      if (deodmh.gt.zero) p2co(nodm1,nodm2)=deodmh
      if (deodmh.gt.zero) p2co(nodm2,nodm1)=deodmh
      if (godmh.gt.zero) p3co(nodm1,nodm2)=godmh
      if (godmh.gt.zero) p3co(nodm2,nodm1)=godmh
      end do
**********************************************************************
*                                                                    *
*     Read in valency angle and conjugation type data                *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990)nvaty
      iline=iline+1
      if (nvaty.gt.nvatym)
     $stop 'Maximum nr. of valency angle types exceeded'
      do i1=1,nvaty
      read (4,1500,end=990,err=990)nvs(i1,1),nvs(i1,2),
     $nvs(i1,3),th0(i1),vka(i1),vka3(i1),vka8(i1),vkac(i1),vkap(i1),
     $vval2(i1)
      iline=iline+1
      end do
**********************************************************************
*                                                                    *
*     Read in torsion angle type data                                *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990)ntoty
      iline=iline+1
      if (ntoty.gt.ntotym)
     $stop 'Maximum nr. of torsion angle types exceeded'
      do i1=1,ntoty
      read (4,1600,end=990,err=990)nts(i1,1),nts(i1,2),nts(i1,3),
     $nts(i1,4),v1(i1),
     $v2(i1),v3(i1),v4(i1),vconj(i1),v2bo(i1),v3bo(i1)
      iline=iline+1
      end do
**********************************************************************
*                                                                    *
*     Read in hydrogen bond type data                                *
*                                                                    *
**********************************************************************
      read (4,1100,end=990,err=990)nhbty
      iline=iline+1
      if (nhbty.gt.nhbtym)
     $stop 'Maximum nr. of hydrogen bond types exceeded'
      do i1=1,nhbty
      read (4,1500,end=990,err=990)nhbs(i1,1),nhbs(i1,2),
     $nhbs(i1,3),rhb(i1),dehb(i1),vhb1(i1),vhb2(i1)
      iline=iline+1
      end do
**********************************************************************
*                                                                    *
*     Calculate vdWaals interaction parameters                       *
*                                                                    *
**********************************************************************
      do i1=1,nso
      do i2=1,nso
      rr=(rvdw(i1)+rvdw(i2))
      rr2=rr*rr
      eps2=sqrt(eps(i1)*eps(i2))
      rr6=rr2*rr2*rr2
      pvdw1(i1,i2)=eps2*rr6*rr6
      pvdw1(i2,i1)=eps2*rr6*rr6
      pvdw2(i1,i2)=2.0*eps2*rr6
      pvdw2(i2,i1)=2.0*eps2*rr6
      end do
      end do
**********************************************************************
*                                                                    *
*     Error part                                                     *
*                                                                    *
**********************************************************************
      goto 999
  990 write (*,*)'Error or end-of-file reading unit 4 on line:',iline
      stop
  999 continue
      close(4)
**********************************************************************
*                                                                    *
*     Format part                                                    *
*                                                                    *
**********************************************************************
 1100 format (i3,2x,a2,3x,3d22.15)
 1200 format (1x,a2,10f9.4)
 1250 format (3x,10f9.4)
 1300 format (f10.4)
 1400 format (2i3,8f9.4)
 1450 format (6x,8f9.4)
 1500 format (3i3,7f9.4)
 1600 format (4i3,7f9.4)
      return
      end
**********************************************************************
***********************************************************************

      subroutine mdsav(node)

***********************************************************************
#include "cbka.blk"
#include "cbkabo.blk"
#include "cbkatomcoord.blk"
#include "cbkbo.blk"
#include "cbkc.blk"
#include "cbkch.blk"
#include "cbkconst.blk"
#include "cbkdistan.blk"
#include "cbkenergies.blk"
#include "cbkia.blk"
#include "cbkinit.blk"
#include "cbklonpar.blk"
#include "cbknubon2.blk"
#include "cbkqa.blk"
#include "cbktregime.blk"
#include "cbksrtbon1.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"

      dimension idum(mbond+3),bodum(mbond+3),qat2(2)
      character*25 qfileh
      character*33 qfile2
      character*4 qext
      character*6 qmdfi
      character *7 var
      character *3 qat2,pepname
      character *1 qrtemp
************************************************************************
*                                                                      *
*     Save coordinates, velocities and accelerations of MD-system      *
*                                                                      *
************************************************************************
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In mdsav'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

************************************************************************
c
c     This is just for test purposes
c
************************************************************************
c$$$      write(6,*) '***************************'
c$$$      write(6,*) 'mdsav node number is ',node
c$$$      write(6,*) '***************************'
      return


      qfileh='Unknown'
      qmdfi='moldyn'
      pepname='   '
      ipeptide=0
      if (ni.eq.2) qmdfi='molsav'

      if (iopt.eq.0) then

      do i1=1,mbond+3
      idum(i1)=nzero
      bodum(i1)=zero
      end do
C      if (napp.eq.1) 
C     $open (7,file='fort.7',status='unknown',access='append')
      if (napp.ne.1) 
     $open (7,file='fort.7',status='unknown')
      nsbmaxh=5*((nsbmax/5)+1)
      write (7,100)na,qmol,mdstep,nsbmaxh
      if (nbiolab.eq.1) write (67,101)na,qmol
      do i1=1,na
      bosum=0.0
      do i3=1,nsbmax
      if (iag(i1,2+i3).gt.0) bosum=bosum+bo(nubon1(i1,i3))
      end do
      if (nsbmax.lt.5) then
      write (7,200)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2)),
     $(idum(i2),i2=1,5-iag(i1,2)),
     $iag(i1,3+mbond),(bo(nubon1(i1,i2)),i2=1,iag(i1,2)),
     $(bodum(i2),i2=1,5-iag(i1,2)),abo(i1),vlp(i1),ch(i1)
      if (nbiolab.eq.1) then      !Delphi-connection table output
      write (67,201)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2))
      end if
      else if (nsbmax.lt.10) then
      write (7,210)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2)),
     $(idum(i2),i2=1,10-iag(i1,2)),
     $iag(i1,3+mbond),(bo(nubon1(i1,i2)),i2=1,iag(i1,2)),
     $(bodum(i2),i2=1,10-iag(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbmax.lt.15) then
      write (7,220)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2)),
     $(idum(i2),i2=1,15-iag(i1,2)),
     $iag(i1,3+mbond),(bo(nubon1(i1,i2)),i2=1,iag(i1,2)),
     $(bodum(i2),i2=1,15-iag(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbmax.lt.20) then
      write (7,230)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2)),
     $(idum(i2),i2=1,20-iag(i1,2)),
     $iag(i1,3+mbond),(bo(nubon1(i1,i2)),i2=1,iag(i1,2)),
     $(bodum(i2),i2=1,20-iag(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbmax.lt.25) then
      write (7,240)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2)),
     $(idum(i2),i2=1,25-iag(i1,2)),
     $iag(i1,3+mbond),(bo(nubon1(i1,i2)),i2=1,iag(i1,2)),
     $(bodum(i2),i2=1,25-iag(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbmax.gt.25) then
      write (7,250)i1,iag(i1,1),(iag(i1,2+i2),i2=1,iag(i1,2)),
     $(idum(i2),i2=1,35-iag(i1,2)),
     $iag(i1,3+mbond),(bo(nubon1(i1,i2)),i2=1,iag(i1,2)),
     $(bodum(i2),i2=1,35-iag(i1,2)),abo(i1),vlp(i1),ch(i1)
      end if
      end do
      boss=zero
      vlps=0.0
C      if (napp.eq.1)
C     $open (8,file='fort.8',status='unknown',access='append')
      if (napp.ne.1)
     $open (8,file='fort.8',status='unknown')
      nsbmaxh=5*((nsbma2/5)+1)
      write (8,100)na,qmol,mdstep,nsbmaxh
      chsum=0.0
      do i1=1,na
      bosum=0.0
      do i3=1,nsbma2
      if (ia(i1,2+i3).gt.0) bosum=bosum+bo(nubon2(i1,i3))
      end do
      if (nsbma2.lt.5) then
      write (8,200)i1,ia(i1,1),(ia(i1,2+i2),i2=1,ia(i1,2)),
     $(idum(i2),i2=1,5-ia(i1,2)),
     $ia(i1,3+mbond),(bo(nubon2(i1,i2)),i2=1,ia(i1,2)),
     $(bodum(i2),i2=1,5-ia(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbma2.lt.10) then
      write (8,210)i1,ia(i1,1),(ia(i1,2+i2),i2=1,ia(i1,2)),
     $(idum(i2),i2=1,10-ia(i1,2)),
     $ia(i1,3+mbond),(bo(nubon2(i1,i2)),i2=1,ia(i1,2)),
     $(bodum(i2),i2=1,10-ia(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbma2.lt.15) then
      write (8,220)i1,ia(i1,1),(ia(i1,2+i2),i2=1,ia(i1,2)),
     $(idum(i2),i2=1,15-ia(i1,2)),
     $ia(i1,3+mbond),(bo(nubon2(i1,i2)),i2=1,ia(i1,2)),
     $(bodum(i2),i2=1,15-ia(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbma2.lt.20) then
      write (8,230)i1,ia(i1,1),(ia(i1,2+i2),i2=1,ia(i1,2)),
     $(idum(i2),i2=1,20-ia(i1,2)),
     $ia(i1,3+mbond),(bo(nubon2(i1,i2)),i2=1,ia(i1,2)),
     $(bodum(i2),i2=1,20-ia(i1,2)),abo(i1),vlp(i1),ch(i1)
      else if (nsbma2.lt.25) then
      write (8,240)i1,ia(i1,1),(ia(i1,2+i2),i2=1,ia(i1,2)),
     $(idum(i2),i2=1,25-ia(i1,2)),
     $ia(i1,3+mbond),(bo(nubon2(i1,i2)),i2=1,ia(i1,2)),
     $(bodum(i2),i2=1,25-ia(i1,2)),abo(i1),vlp(i1),ch(i1)
      end if
      boss=boss+bosum/2.0
      vlps=vlps+vlp(i1)
      chsum=chsum+ch(i1)
      end do
      write (7,*)2.0*boss,vlps,2.0*boss+2.0*vlps,chsum
      close(8)
      close(7)

      end if

      if (noutpt.eq.0) then
      write (var,'(f7.4)')float(mdstep/nsav)/1d4
      if (ni.eq.0) open (unit=67,file=qmdfi//var(3:7),
     $status='unknown')
      write (67,300)qmol
      do i1=1,na
      write (67,400)i1,qa(i1),(c(i1,i2),i2=1,3)
      end do
      write (67,*)
      close(67)
      end if

      if (noutpt.eq.2) then
C      open (88,file='moldyn.bgf',status='unknown',access='append')
      call writebgf(88)
      close (88)
      end if

      if ((ni.eq.1.and.iopt.eq.0).or.(ni.eq.1.and.iopt.eq.1.and.
     $iflga.eq.1)) then
      qrtemp=qr
      if (qr.eq.'I') qr='C'
      if (qfileh.eq.' ') then
      write (*,*)'Warning: no file name given; use Unknown' 
      qfileh='Unknown'
      end if
      qfile2=qfileh
      if (imodfile.eq.0) then
      istart=1
      qstrana1(1:25)=qfileh
      call stranal(istart,iend,vout,iout,1)
      qfile2=qfileh(istart:iend-1)//".geo"
      end if
      call writegeo(98)

      if (imodfile.eq.1.or.iopt.eq.0) then
      open (88,file=qfile2,status='unknown')
      call writegeo(88)
      close (88)
      end if 

      qr=qrtemp

      if (iopt.eq.0) then

      do i1=1,na
      write (56,410) i1,ch(i1)
      write (55,410) i1,chgbgf(i1)
      end do
**********************************************************************
*                                                                    *
*     Write .pdb output file                                         *
*                                                                    *
**********************************************************************
      open (unit=47,file='output.pdb',status='unknown')
      do i1=1,na
      write (47,412)'ATOM  ',i1,qa(i1),pepname,ipeptide,c(i1,1),
     $c(i1,2),c(i1,3),1.0,2.2,qa(i1)
      end do
      write (47,*) 'TER'
      write (47,*) 'END'
      close (47)

      if (nsurp.eq.0) then
      if (kx.gt.0.or.ky.gt.0.or.kz.gt.0) then
      qrtemp=qr
**********************************************************************
*                                                                    *
*     Write crystal structure including periodic images              *
*                                                                    *
**********************************************************************
*     mux=(1+kx+kx)
*     muy=(1+ky+ky)
*     muz=(1+kz+kz)
*     qr='F'
*     write (86,'(2x,a1,1x,a60)')qr,qmol
*     qr=qrtemp
*     write (86,'(3f10.4)')mux*axiss(1),muy*axiss(2),muz*axiss(3)
*     write (86,'(3f10.4)')angle(1),angle(2),angle(3)
*     do i1=1,na
*     write (86,'(i4,1x,a2,3x,3d22.15)')i1,qa(i1),(c(i1,i2),i2=1,3)
*     end do
*     nhulp=na+1
*     do k1=-kx,kx
*     do k2=-ky,ky
*     do k3=-kz,kz
*     if (k1.ne.0.or.k2.ne.0.or.k3.ne.0) then
*     do i1=1,na
*     cx=c(i1,1)+k1*tm11
*     cy=c(i1,2)+k1*tm21+k2*tm22
*     cz=c(i1,3)+k1*tm31+k2*tm32+k3*tm33
*     write (86,'(i4,1x,a2,3x,3d22.15)')nhulp,qa(i1),cx,cy,cz
*     nhulp=nhulp+1
*     end do
*     end if
*     end do
*     end do
*     end do
*     write (86,*)
**********************************************************************
*                                                                    *
*     Write crystal structure with extra unit cells                  *
*                                                                    *
**********************************************************************
      mux=1+iexx
      muy=1+iexy
      muz=1+iexz
      qr='F'
      write (85,'(2x,a1,1x,a60)')qr,qmol
      qr=qrtemp
      write (85,'(3f10.4)')mux*axiss(1),muy*axiss(2),muz*axiss(3)
      write (85,'(3f10.4)')angle(1),angle(2),angle(3)
      do i1=1,na
      write (85,'(i4,1x,a2,3x,3d22.15)')i1,qa(i1),(c(i1,i2),i2=1,3)
      end do
      nhulp=na+1
      do k1=0,iexx
      do k2=0,iexy
      do k3=0,iexz
      if (k1.ne.0.or.k2.ne.0.or.k3.ne.0) then
      do i1=1,na
      cx=c(i1,1)+k1*tm11
      cy=c(i1,2)+k1*tm21+k2*tm22
      cz=c(i1,3)+k1*tm31+k2*tm32+k3*tm33
      write (85,'(i4,1x,a2,3x,3d22.15)')nhulp,qa(i1),cx,cy,cz
      nhulp=nhulp+1
      end do
      end if
      end do
      end do
      end do
      write (85,*)

      end if
      end if
      end if

      end if

      if (ni.eq.0.or.ni.eq.2) then
********************************************************************** 
*                                                                    *
*     Write ASCII trajectory file                                    *
*                                                                    *
********************************************************************** 
      if (ni.eq.0) open(unit=66,file=qmdfi//'.vel',status='unknown')
      if (ni.eq.2) then
      write (var,'(f7.4)')float(mdstep/nsav3)/1d4
      open (unit=66,file=qmdfi//var(3:7),status='unknown')
      end if
      write (66,500)axis(1),axis(2),axis(3)
      write (66,550)angle(1),angle(2),angle(3)
      write (66,600)na,((c(i,j),j=1,3),qlabel(i),i=1,na)
      write (66,700)((vel(j,i),j=1,3),i=1,na)
      write (66,800)((accel(j,i),j=1,3),i=1,na)
      write (66,900)((aold(j,i),j=1,3),i=1,na)
      write (66,1000)tempmd
      write (66,1050)
      close (66)
      end if
      if (ni.ne.2.and.iopt.eq.0) then

C      open (unit=68,file='xmolout',status='unknown',access='append')
      write (68,1200)na
      write (68,1300)qmol,mdstep+nit+nprevrun,estrc,
     $axis(1),axis(2),axis(3),angle(1),angle(2),angle(3)
      do i1=1,na
      if (ixmolo.eq.0) write (68,1400)qa(i1),(c(i1,i2),i2=1,3)
      if (ixmolo.eq.1) write (68,1400)qa(i1),(c(i1,i2),i2=1,3),
     $(vel(i2,i1)/1e+10,i2=1,3),iag(i1,3+mbond)
      if (ixmolo.eq.2) write (68,1401)qa(i1),(c(i1,i2),i2=1,3),
     $iag(i1,3+mbond)
      end do
      close (68)

      if (itrout.ne.0) then
C      open (unit=69,file='xmolout2',status='unknown',access='append')
      write (69,1200)na
      write (69,1300)qmol,mdstep+nit+nprevrun,estrc,
     $axis(1),axis(2),axis(3),angle(1),angle(2),angle(3)
      do i1=1,na
      if (ixmolo.eq.0) write (69,1400)qa(i1),(cp(i1,i2),i2=1,3)
      if (ixmolo.eq.1) write (69,1400)qa(i1),(cp(i1,i2),i2=1,3),
     $(vel(i2,i1)/1e+10,i2=1,3),iag(i1,3+mbond)
      if (ixmolo.eq.2) write (68,1401)qa(i1),(c(i1,i2),i2=1,3),
     $iag(i1,3+mbond)
      end do
      close (69)
      end if

      call molanal
      end if
********************************************************************** 
*                                                                    *
*     Generate BIOGRAF output-file                                   *
*                                                                    *
********************************************************************** 
      if ((ni.eq.1.and.iopt.eq.0).or.(ni.eq.1.and.iopt.eq.1.and.
     $iflga.eq.1)) then

      if (qfileh.eq.' ') then
      write (*,*)'Warning: no file name given; use Unknown' 
      qfileh='Unknown'
      end if
      qfile2=qfileh
      if (imodfile.eq.0) then
      istart=1
      qstrana1(1:25)=qfileh
      call stranal(istart,iend,vout,iout,1)
      qfile2=qfileh(istart:iend-1)//".bgf"
      end if
      call writebgf(90)

      if (imodfile.eq.1.or.iopt.eq.0) then
      open (88,file=qfile2,status='unknown')
      call writebgf(88)
      close (88)
      end if

      end if
     
      return
********************************************************************** 
*                                                                    *
*     Format part                                                    *
*                                                                    *
********************************************************************** 
  100 format (i4,1x,a40,'Iteration:',i8,' #Bonds:',i4)       
  101 format (i3,2x,a40)
  200 format (8i4,8f7.3)
  201 format (8i3)
  210 format (13i4,13f7.3)
  220 format (18i4,18f7.3)
  230 format (23i4,23f7.3)
  240 format (28i4,28f7.3)
  250 format (38i4,38f7.3)
  300 format (2x,a1,1x,a60)
  301 format (2x,a1,1x,f6.2,a60)
  302 format (2x,a1,1x,2f6.2,a60)
  310 format (2x,a1,1x,a60)
  320 format (3f10.4)
  400 format (i4,1x,a2,3x,3(d21.14,1x),1x,a5,1x,i5)
  410 format (i4,f12.6)
  412 format(A6,I5,1x,A2,3x,A3,2x,i4,4x,3f8.3,f6.2,f6.2,4x,2x,A6)
  500 format (1x,'Lattice parameters:',/(3f15.8))
  550 format (3f15.8)
  600 format (i4,1x,'Atom coordinates (Angstrom):',/
     $(3d24.15,1x,a5))
  700 format (1x,'Atom velocities (Angstrom/s):',/(3d24.15))
  800 format (1x,'Atom accelerations (Angstrom/s**2):',/(3d24.15))
  900 format (1x,'Previous atom accelerations:',/(3d24.15))
 1000 format (1x,'MD-temperature (K):',/(1d24.15))
 1050 format (1x,'Connections, bond orders and lone pairs:')
 1100 format (8i3,8f8.4)
 1200 format (i4)
 1300 format (a40,i6,f12.4,6f7.2)
 1400 format (a2,3f10.5,3f15.5,i6)
 1401 format (a2,3f10.5,i6)
 1500 format ('BIOGRF',i4)
 1600 format ('XTLGRF',i4)
 1700 format ('DESCRP ',a60)
 1800 format ('REMARK ',a60)
 1900 format ('FFIELD ',a40)
 2000 format ('RUTYPE ',a40)
 2100 format ('CRYSTX ',6f11.5)
 2200 format ('CELLS ',6i5)
 2300 format ('#              At1 At2   R12    Force1  Force2  ',
     $'dR12/dIteration(MD only)')
 2400 format ('BOND RESTRAINT ',2i4,f8.4,f8.2,f8.5,f10.7)
 2500 format ('#               At1 At2 At3 Angle   Force1  Force2',
     $'  dAngle/dIteration (MD only)')
 2600 format ('ANGLE RESTRAINT ',3i4,2f8.2,f8.4,f9.6)
 2700 format ('#                 At1 At2 At3 At3 Angle   Force1  ',
     $'Force2  dAngle/dIteration (MD only)')
 2800 format ('TORSION RESTRAINT ',4i4,2f8.2,f8.4,f9.6)
 2900 format ('FORMAT ATOM   (a6,1x,i5,1x,a5,1x,a3,1x,a1,1x,a5,',
     $'3f10.5,1x,a5,i3,i2,1x,f8.5)')
 3000 format ('HETATM',1x,i5,1x,a5,1x,a3,1x,a1,1x,a5,3f10.5,1x,
     $a5,i3,i2,1x,f8.5)
 3100 format ('FORMAT CONECT (a6,12i6)')
 3200 format ('CONECT',12i6)
 3300 format ('UNIT ENERGY   kcal')
 3400 format ('ENERGY',5x,f14.6)
 3500 format ('END')
      end

************************************************************************
************************************************************************

      subroutine readc

************************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkcha.blk"
#include "cbkconst.blk"
#include "cbkdistan.blk"
#include "cbkinit.blk"
#include "cbktregime.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"

      character*6 qident
      character*20 qhulp
*     dimension qident(100)
************************************************************************
*                                                                      *
*     Read control file                                                *
*                                                                      *
************************************************************************
c$$$c      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$c      write (65,*) 'In readc'
c$$$c      call timer(65)
c$$$c      close (65)
c$$$c      end if
      if (mdstep.gt.0.or.nit.gt.0) nmmsav=nmm
************************************************************************
*                                                                      *
*     Set default values                                               *
*                                                                      *
************************************************************************
      nreac=0
      axis1=200.0d0
      axis2=200.0d0
      axis3=200.0d0
      cutof2=0.001d0
      cutof3=0.300d0
      tsetor=298.0d0
      tset2=298.0d0
      pset=0.0d0
      tincr=0.0d0
      tstep=0.5d0
*     swa=0.0   !Moved to force field
*     swb=12.5  !Moved to force field
      taut=2.5d0
      taut2=2.5d0
      ndtau=50000
      taup=500.0d0
      vqnd=100.0d0
      errnh=1.0d0
      range=2.5d0
      maxstp=1000
      nequi=0
      nmethod=3
      ncha=3
      ncha2=1
      nchaud=1
      nvlist=25
      nrep1=5
      nsav=50
      icheck=0
      ivels=0
      itfix=0
      ncontrol=25
      noutpt=0
      napp=0
      nsurp=0
      ncons=2
      nrand=0
      nmm=0
      endpo=1.0d0
      endpo2=1.0d0
      nfc=50
      nsav2=50
      nmmax=50
      i5758=0
      parc1=1.0d0
      parc2=0.001d0
      icell=0
      ingeo=1
      ccpar=1.0005d0
      icelo2=0
      nrdd=0
      nrddf=200000
      nbiolab=0
c      ngeofor=0
      nincrop=0
      accerr=2.50d0
      vrange=2.50d0
      vlbora=5.00d0
      nsav3=1000
      nhop2=25
      nprevrun=0
      ndebug=0
      volcha=10.00d0
      ixmolo=0
      inpt=0
      iconne=0
      imolde=0
      ianaly=0
      icentr=0
      itrans=0
      itrout=0
      tpnrad=300.0d0
      ityrad=3
      iexx=1
      iexy=1
      iexz=1
      syscha=0.00d0
      inmov1=0
      inmov2=0
      vfield=0.00d0
      itstep=0
      ifreq=0
      isymm=1
      icpres=0
      delvib=0.0001d0
c     shock variables
      shock_vel = 2.d0 ! impact velocity for shock simulations (nm/ps)
      shock_z_sep = 10.0d0 ! separation z value to apply initial velocities in shocks
      ishock_type = 0.0d0 ! shock type. 0: simple impact; 1: compressing c axis
c     Hugoniostat variables
      Hug_E0 = 0.d0 ! Reference energy
      Hug_P0 = 0.d0 ! Reference pressure
      Hug_V0 = 0.d0 ! Reference volume
c     Shear flow simulations for viscosity
      xImpVcm = 1.d0 ! velocity applied in shear simulations (in nm/ps), left half mover at -xImpVcm and right at +xImpVcm
c$$$************************************************************************
c$$$*                                                                      *
c$$$*     Read control-file                                                *
c$$$*                                                                      *
c$$$************************************************************************
c$$$      open (10,file='control',status='old')
c$$$   10 read (10,'(a20)',end=20,err=30)qhulp
c$$$      if (qhulp(1:1).eq.'#') goto 10
c$$$      read (qhulp,*,err=30)vhulp
c$$$      read (qhulp,'(8x,a6)',err=30)qident
c$$$      if (qident.eq.'Hug_V0') Hug_P0=vhulp
c$$$      if (qident.eq.'Hug_P0') Hug_V0=vhulp
c$$$      if (qident.eq.'Hug_E0') Hug_E0=vhulp
c$$$      if (qident.eq.'shea_v') xImpVcm=vhulp
c$$$      if (qident.eq.'shok_t') ishock_type=int(vhulp)
c$$$      if (qident.eq.'shok_z') shock_z_sep=vhulp
c$$$      if (qident.eq.'shok_v') shock_vel=vhulp
c$$$      if (qident.eq.'nreac') nreac=int(vhulp)
c$$$      if (qident.eq.'axis1') axis1=vhulp
c$$$      if (qident.eq.'axis2') axis2=vhulp
c$$$      if (qident.eq.'axis3') axis3=vhulp
c$$$      if (qident.eq.'cutof2') cutof2=vhulp
c$$$      if (qident.eq.'cutof3') cutof3=vhulp
c$$$      if (qident.eq.'mdtemp') tsetor=vhulp
c$$$      if (qident.eq.'mdtem2') tset2=vhulp
c$$$      if (qident.eq.'mdpres') pset=vhulp*0.001
c$$$      if (qident.eq.'tincr') tincr=vhulp
c$$$      if (qident.eq.'tstep') tstep=vhulp
c$$$*     if (qident.eq.'lowtap') swa=vhulp !Moved to force field
c$$$*     if (qident.eq.'uptap') swb=vhulp  !Moved to force field
c$$$      if (qident.eq.'tdamp1') taut=vhulp
c$$$      if (qident.eq.'tdamp2') taut2=vhulp
c$$$      if (qident.eq.'ntdamp') ndtau=int(vhulp)
c$$$      if (qident.eq.'pdamp1') taup=vhulp
c$$$      if (qident.eq.'tdhoov') vqnd=vhulp
c$$$      if (qident.eq.'achoov') errnh=vhulp/100.0
c$$$      if (qident.eq.'range') range=vhulp
c$$$      if (qident.eq.'nmdit') maxstp=int(vhulp)
c$$$      if (qident.eq.'nmdeqi') nequi=int(vhulp)
c$$$      if (qident.eq.'imdmet') nmethod=int(vhulp)
c$$$      if (qident.eq.'icharg') ncha=int(vhulp)
          nchaold=ncha
c$$$      if (qident.eq.'ichaen') ncha2=int(vhulp)
c$$$      if (qident.eq.'ichupd') nchaud=int(vhulp)
c$$$      if (qident.eq.'iout1') nrep1=int(vhulp)
c$$$      if (qident.eq.'iout2') nsav=int(vhulp)
c$$$      if (qident.eq.'icheck') ntest=int(vhulp)
c$$$      if (qident.eq.'ivels') nvel=int(vhulp)
c$$$      if (qident.eq.'itfix') ntscale=int(vhulp)
c$$$      if (qident.eq.'irecon') ncontrol=int(vhulp)
c$$$      if (qident.eq.'iout3') noutpt=int(vhulp)
c$$$      if (qident.eq.'iappen') napp=int(vhulp)
c$$$      if (qident.eq.'isurpr') nsurp=int(vhulp)
c$$$      if (qident.eq.'itdmet') ncons=int(vhulp)
c$$$      if (qident.eq.'iravel') nrand=int(vhulp)
c$$$      if (qident.eq.'imetho') nmm=int(vhulp)
c$$$      if (qident.eq.'endmm') endpo=vhulp
          endpoold=endpo
c$$$      if (qident.eq.'endmd') endpo2=vhulp
c$$$      if (qident.eq.'imaxmo') nfc=int(vhulp)
          nfcold=nfc
c$$$      if (qident.eq.'iout4') nsav2=int(vhulp)
c$$$      if (qident.eq.'imaxit') nmmax=int(vhulp)
          nmmaxold=nmmax
c$$$      if (qident.eq.'iout5') i5758=int(vhulp)
c$$$      if (qident.eq.'parsca') parc1=vhulp
c$$$      if (qident.eq.'parext') parc2=vhulp
c$$$      if (qident.eq.'icelop') icell=int(vhulp)
          icellold=icell
c$$$      if (qident.eq.'igeopt') ingeo=int(vhulp)
c$$$      if (qident.eq.'celopt') ccpar=vhulp
c$$$      if (qident.eq.'icelo2') icelo2=int(vhulp)
          icelo2old=icelo2
c$$$      if (qident.eq.'ideve1') nrdd=int(vhulp)
c$$$      if (qident.eq.'ideve2') nrddf=int(vhulp)
c$$$      if (qident.eq.'ibiola') nbiolab=int(vhulp)
c$$$c      if (qident.eq.'igeofo') ngeofor=int(vhulp)
c$$$      if (qident.eq.'iincop') nincrop=int(vhulp)
c$$$      if (qident.eq.'accerr') accincr=vhulp
c$$$      if (qident.eq.'iout6') nsav3=int(vhulp)
c$$$      if (qident.eq.'irten') nhop2=int(vhulp)
c$$$      if (qident.eq.'npreit') nprevrun=int(vhulp)
c$$$      if (qident.eq.'idebug') ndebug=int(vhulp)
c$$$      if (qident.eq.'volcha') volcha=vhulp
c$$$      if (qident.eq.'ixmolo') ixmolo=int(vhulp)
c$$$      if (qident.eq.'inpt') inpt=int(vhulp)
c$$$      if (qident.eq.'iconne') iconne=int(vhulp)
c$$$      if (qident.eq.'imolde') imolde=int(vhulp)
c$$$      if (qident.eq.'ianaly') ianaly=int(vhulp)
c$$$      if (qident.eq.'icentr') icentr=int(vhulp)
c$$$      if (qident.eq.'itrans') itrans=int(vhulp)
c$$$      if (qident.eq.'itrout') itrout=int(vhulp)
c$$$      if (qident.eq.'nvlist') nvlist=int(vhulp)
c$$$      if (qident.eq.'vrange') vrange=vhulp
c$$$      if (qident.eq.'vlbora') vlbora=vhulp
c$$$      if (qident.eq.'tpnrad') tpnrad=vhulp
c$$$      if (qident.eq.'ityrad') ityrad=int(vhulp)
c$$$      if (qident.eq.'iexx') iexx=int(vhulp)
c$$$      if (qident.eq.'iexy') iexy=int(vhulp)
c$$$      if (qident.eq.'iexz') iexz=int(vhulp)
c$$$      if (qident.eq.'syscha') syscha=vhulp
c$$$      if (qident.eq.'inmov1') inmov1=int(vhulp)
c$$$      if (qident.eq.'inmov2') inmov2=int(vhulp)
c$$$      if (qident.eq.'itstep') itstep=int(vhulp)
c$$$      if (qident.eq.'ifreq') ifreq=int(vhulp)
c$$$      if (qident.eq.'isymm') isymm=int(vhulp)
c$$$      if (qident.eq.'icpres') icpres=int(vhulp)
c$$$      if (qident.eq.'delvib') delvib=vhulp
c$$$      goto 10
c$$$   20 continue
      close (10)
      axis(1)=axis1
      axis(2)=axis2
      axis(3)=axis3
      if (axiss(1).gt.zero) then
      axis(1)=axiss(1)
      axis(2)=axiss(2)
      axis(3)=axiss(3)
      end if
      if (tincr.lt.0.0001.and.tincr.gt.-0.0001) tset=tsetor
      iequi=1
      if (nequi.gt.0) iequi=0
      if (iopt.eq.1.and.napp.eq.1) then
      stop 'No fort.7 and fort.8 append with iopt=1 !'
      end if
      if (mdstep.gt.0.or.nit.gt.0) nmm=nmmsav
      if (mdstep.gt.0.and.itstep.eq.1) then
      tstepmax=tstep
      tstep=tstep*(tsetor/tempmd)
      if (tstep.gt.tstepmax) tstep=tstepmax
      end if
      tstep=1.0d-15*tstep
      taus=taut
      taut=1.0d-15*taut
      taut2=1.0d-15*taut2
      taup=1.0d-15*taup
      ts2=tstep/2.0
      ts22=tstep*ts2
      return
   30 continue
      write (*,*)'Error reading control-file'
      stop 'Error reading control-file'
************************************************************************
*                                                                      *
*     Format part                                                      *
*                                                                      *
************************************************************************
 1050 format (f7.3)
 1055 format (f7.4)
 1056 format (f9.4)
 1060 format (i8)
 1070 format (f7.5)
      end
************************************************************************
************************************************************************

      subroutine staint

************************************************************************
#include "cbka.blk"
#include "cbkdcell.blk"
#include "cbkqa.blk"
#include "control.blk"
#include "small.blk"
#include "cbkc.blk"
#include "cbkconst.blk"
      dimension bvt(nat,4)
************************************************************************
*                                                                      *
*     Generate cartesian coordinates from internal coordinate input    *
*                                                                      *
************************************************************************
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In staint'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      k=0
   10 read (3,1200,end=20,err=20)(ijk(k+1,k1),k1=1,3),k2,qa(k+1),
     $bvt(k+1,3),bvt(k+1,2),bvt(k+1,1)
      qlabel(k+1)=qa(k+1)
      qresi1(k+1)='   '
      qresi2(k+1)=' '
      qresi3(k+1)='     '
      qffty(k+1)='     '
      if (k2.ne.k+1) then
      write (*,*)'Wrong order in internal coordinates at atom:',k2
      goto 20
*     stop 'Wrong order in internal coordinates'
      end if
      k=k+1
      if (k.gt.nat) then
      write (*,*)na,nat
      stop 'Maximum number of atoms exceeded'
      end if
      goto 10
   20 continue
      na=k
      
************************************************************************
*                                                                      *
*     CALCULATION OF CARTESIAN COORDINATES FROM INTERNAL COORDINAATES  *
*                                                                      *
************************************************************************

   12 C(1,1)=ZERO
      C(1,2)=ZERO
      C(1,3)=ZERO
      C(2,1)=BVT(2,1)
      C(2,2)=ZERO
      C(2,3)=ZERO
      HR=(BVT(3,2)-90.0D0)*DGRRDN
      C(3,1)=C(2,1)+BVT(3,1)*SIN(HR)
      C(3,2)=BVT(3,1)*COS(HR)
      C(3,3)=ZERO
      DO 32 K1=4,NA
      J=IJK(K1,2)
      KB=K1-1
      XH=C(J,1)
      YH=C(J,2)
      ZH=C(J,3)
      DO 13 K2=1,KB
      C(K2,1)=C(K2,1)-XH
      C(K2,2)=C(K2,2)-YH
      C(K2,3)=C(K2,3)-ZH
      DO 13 K3=1,3
   13 IF (ABS(C(K2,K3)).LT.1.0D-15) C(K2,K3)=1.0D-15
      K=IJK(K1,3)
      P2=C(K,2)*C(K,2)+C(K,3)*C(K,3)
      IF (P2.NE.ZERO) THEN
      P=SQRT(P2)
      Q=SQRT(C(K,1)*C(K,1)+P2)
      SA=C(K,2)/P
      CA=C(K,3)/P
      SB=-C(K,1)/Q
      CB=P/Q
      ELSE
      SA=ZERO
      CA=ONE
      SB=ONE
      CB=ZERO
      ENDIF
      DO 16 K2=1,KB
      AZ=C(K2,1)
      BZ=C(K2,2)
      C(K2,1)=AZ*CB+BZ*SB*SA+C(K2,3)*SB*CA
      C(K2,2)=BZ*CA-C(K2,3)*SA
   16 C(K2,3)=-AZ*SB+BZ*CB*SA+C(K2,3)*CB*CA
      IF (C(K,3).LE.ZERO) THEN
      DO 17 K2=1,KB
   17 C(K2,3)=-C(K2,3)
      ENDIF
      I=IJK(K1,1)
      IF (1.0D5*ABS(C(I,1)).LE.ABS(C(I,2))) THEN
      T1=HALF*PI
      ELSE
      YX=ABS(C(I,2)/C(I,1))
      T1=ATAN(YX)
      ENDIF
      IF (C(I,1).GE.ZERO.AND.C(I,2).LT.ZERO) T1=TWO*PI-T1
      IF (C(I,1).LT.ZERO.AND.C(I,2).GE.ZERO) T1=PI-T1
      IF (C(I,1).LT.ZERO.AND.C(I,2).LT.ZERO) T1=T1+PI
      DO 31 K2=1,KB
      IF (C(K2,1).EQ.ZERO.AND.C(K2,2).EQ.ZERO) GOTO 31
      IF (1.0D5*ABS(C(K2,1)).LT.ABS(C(K2,2))) THEN
      T2=HALF*PI
      ELSE
      YX=ABS(C(K2,2)/C(K2,1))
      T2=ATAN(YX)
      ENDIF
      IF (C(K2,1).GE.ZERO.AND.C(K2,2).LT.ZERO) T2=TWO*PI-T2
      IF (C(K2,1).LT.ZERO.AND.C(K2,2).GE.ZERO) T2=PI-T2
      IF (C(K2,1).LT.ZERO.AND.C(K2,2).LT.ZERO) T2=T2+PI
      T3=T2-T1
      IF (T3.LT.ZERO)T3=T3+TWO*PI
      RZ=SQRT(C(K2,1)*C(K2,1)+C(K2,2)*C(K2,2))
      C(K2,1)=RZ*COS(T3)
      C(K2,2)=RZ*SIN(T3)
   31 CONTINUE
      HR=(BVT(K1,2)-90.0D0)*DGRRDN
      HT=BVT(K1,3)*DGRRDN
      CHR=COS(HR)
      C(K1,1)=BVT(K1,1)*CHR*COS(HT)
      C(K1,2)=BVT(K1,1)*CHR*SIN(HT)
   32 C(K1,3)=C(IJK(K1,3),3)+BVT(K1,1)*SIN(HR)
      
      return
 1200 FORMAT(4I3,1X,A2,3F10.5,4X,I1,F10.5)
      end
************************************************************************
************************************************************************

      subroutine outint

************************************************************************
#include "cbka.blk"
#include "cbkabo.blk"
#include "cbkbo.blk"
#include "cbkconst.blk"
#include "cbkia.blk"
#include "cbkinit.blk"
#include "cbknubon2.blk"
#include "cbkqa.blk"
#include "cbktregime.blk"
#include "control.blk"
#include "small.blk"
#include "cbksrtbon1.blk"
************************************************************************
*                                                                      *
*     Output internal coordinates                                      *
*                                                                      *
************************************************************************
      dimension dvdc(3,3),dargdc(3,3)
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In outint'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      write (91,50)qmol
      open (82,file='output.MOP',status='unknown')
      write (82,*)
      write (82,'(a40)')qmol
      write (82,*)
      close (82)

*     IF (NMOLO.GT.1) THEN
*     WRITE(6,*)' OUTPUT INTERNAL COORDINATES NOT POSSIBLE FOR ',
*    $'CALCULATION ON MORE THAN ONE MOLECULE'
*     RETURN
*     END IF

************************************************************************
*                                                                      *
*     Output of internal coordinates.                                  *
*     First 3 atoms of other input file.                               *
*                                                                      *
************************************************************************
      N1=1
      N2=2
      N3=3
C      open (82,file='output.MOP',status='unknown',access='append')
      write(91,100)N1,qa(n1)
      write(82,'(2x,a2,f12.6,i3,f12.6,i3,f12.6,i3,1x,3i4)')qa(n1),
     $zero,nzero,zero,nzero,zero,nzero,nzero,nzero,nzero
      call dista2(n1,n2,rr,dx,dy,dz)
      write(91,200)N1,N2,qa(n2),RR
      write(82,'(2x,a2,f12.6,i3,f12.6,i3,f12.6,i3,1x,3i4)')qa(n2),
     $rr,none,zero,nzero,zero,nzero,n1,nzero,nzero
      close (82)

      call dista2(n2,n3,rr,dx,dy,dz)
      hv=zero
      call calvalres(n1,n2,n3,arg,hv,dvdc,dargdc)
      WRITE(91,300)N1,N2,N3,qa(n3),rdndgr*HV,RR
C      open (82,file='output.MOP',status='unknown',access='append')
      write(82,'(2x,a2,f12.6,i3,f12.6,i3,f12.6,i3,1x,3i4)')qa(n3),
     $rr,none,rdndgr*hv,none,zero,nzero,n2,n1,nzero
      close (82)

      naih=3
      
      do i1=naih+1,na
      bomax=zero
      j1=0
      do i2=1,ia(i1,2)
      iob=ia(i1,2+i2)
      ncubo=nubon2(i1,i2)
      if (bo(ncubo).gt.bomax.and.iob.lt.i1) then
      bomax=bo(ncubo)
      j1=iob
      end if
      end do
      if (j1.eq.0) j1=i1-1
      call dista2(j1,i1,rr,dx,dy,dz)
      
      bomax=zero
      j2=0
      do i2=1,ia(j1,2)
      iob=ia(j2,2+i2)
      ncubo=nubon2(j1,i2)
      if (bo(ncubo).gt.bomax.and.iob.lt.i1.and.
     $abo(iob).gt.bo(ncubo)+0.2) then
      bomax=bo(ncubo)
      j2=iob
      end if
      end do
      if (j2.eq.0) j2=i1-2
      if (j2.eq.j1) j2=j2+1

      call calvalres(j2,j1,i1,arg,hh,dvdc,dargdc)

      bomax=zero
      j3=0
      do i2=1,ia(j2,2)
      iob=ia(j2,2+i2)
      ncubo=nubon2(j2,i2)
      if (bo(ncubo).gt.bomax.and.iob.lt.i1.and.iob.ne.j1) then
      bomax=bo(ncubo)
      j3=iob
      end if
      end do
      if (j3.eq.0) j3=i1-3
      if (j3.eq.j2.and.j3.ne.j1-1) j3=j3+1
      if (j3.eq.j2.and.j3.ne.j1-2) j3=j3+2
      if (j3.eq.j1.and.j3.ne.j2-1) j3=j3+1
      if (j3.eq.j1.and.j3.ne.j2-2) j3=j3+2

      call caltor(j3,j2,j1,i1,ht)

      write(91,400)j3,j2,j1,i1,qa(i1),ht,rdndgr*hh,rr
C      open (82,file='output.MOP',status='unknown',access='append')
      write(82,'(2x,a2,f12.6,i3,f12.6,i3,f12.6,i3,1x,3i4)')qa(i1),
     $rr,none,rdndgr*hh,none,ht,none,j1,j2,j3
      close (82)
      end do

      close(82)
      return
   50 format ('  I',2x,a60)
  100 FORMAT(9X,I3,1x,a2)
  200 FORMAT(6X,2I3,1x,a2,20X,F10.5)
  300 FORMAT(3X,3I3,1x,a2,10X,2F10.5)
  400 FORMAT(4I3,1x,a2,3F10.5)

      end
************************************************************************
************************************************************************

      subroutine outres

************************************************************************
#include "cbka.blk"
#include "cbkbo.blk"
#include "cbkch.blk"
#include "cbkd.blk"
#include "cbkenergies.blk"
#include "cbkh.blk"
#include "cbkimove.blk"
#include "cbkrbo.blk"
#include "cbktorang.blk"
#include "cbktorsion.blk"
#include "cbktregime.blk"
#include "cbkvalence.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
#include "cbkinit.blk"

************************************************************************
*                                                                      *
*     Output molecular data                                            *
*                                                                      *
************************************************************************
      dimension isort(100),iad1(100),iad2(100),iad3(100),iad4(100)
      character*60 qm2
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In outres'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      read (9,100,end=50)idata,qm2
*     if (qm2.ne.qmol) then
*     write (*,*)'Wrong molecule in outres-file'
*     write (*,*)qmol
*     write (*,*)qm2
*     return
*     end if
      do 25 i1=1,idata
      read (9,200)isort(i1),iad1(i1),iad2(i1),iad3(i1),iad4(i1)
      ndata2=ndata2+1

      if (isort(i1).eq.1) then
*     do i2=1,nbon
*     if (ib(i2,2).eq.iad1(i1).and.ib(i2,3).eq.iad2(i1)) then
*     if (iopt.ne.1) write (81,*)iad1(i1),iad2(i1),rbo(i2)
*     caldat(ndata2)=rbo(i2)
*     end if
*     end do
      call dista2(iad1(i1),iad2(i1),dish,dx,dy,dz)
      write (81,*)iad1(i1),iad2(i1),dish
      caldat(ndata2)=dish
      end if

      if (isort(i1).eq.2) then
      do i2=1,nval
      if (iv(i2,2).eq.iad1(i1).and.iv(i2,3).eq.iad2(i1).and.
     $iv(i2,4).eq.iad3(i1)) then
      if (iopt.ne.1) write (81,*)iad1(i1),iad2(i1),
     $iad3(i1),h(i2)*rdndgr
      caldat(ndata2)=h(i2)*rdndgr
      end if
      end do
      end if

      if (isort(i1).eq.3) then
      do i2=1,ntor
      if (it(i2,2).eq.iad1(i1).and.it(i2,3).eq.iad2(i1).and.
     $it(i2,4).eq.iad3(i1).and.it(i2,5).eq.iad4(i1)) then
      if (iopt.ne.1) write (81,*)iad1(i1),iad2(i1),iad3(i1),iad4(i1),
     $abs(thg(i2))
      caldat(ndata2)=abs(thg(i2))
      end if
      end do
      end if

      if (isort(i1).eq.4) then
      if (iopt.ne.1) write (81,*)estrmin
      caldat(ndata2)=estrmin
      end if

      if (isort(i1).eq.5) then
      if (iopt.ne.1) write (81,*)estrmin
      caldat(ndata2)=estrmin
      end if

      if (isort(i1).eq.6) then
      if (iopt.ne.1) write (81,*)iad1(i1),axiss(iad1(i1))
      caldat(ndata2)=axiss(iad1(i1))
      end if

      if (isort(i1).eq.7) then
      if (iopt.ne.1) write (81,*)eco
      caldat(ndata2)=eco
      end if

      if (isort(i1).eq.8) then
      do i2=1,nbon
      if (ib(i2,2).eq.iad1(i1).and.ib(i2,3).eq.iad2(i1)) then
      if (iopt.ne.1) write (81,*)iad1(i1),iad2(i1),bo(i2)
      caldat(ndata2)=bo(i2)
      end if
      end do
      end if

      if (isort(i1).eq.9) then
      if (iopt.ne.1) write (81,*)ch(iad1(i1))
      caldat(ndata2)=ch(iad1(i1))
      end if
      
      if (isort(i1).eq.10) then
      rmsg=0.0
      nmovh=0
      do i2=1,na
      do i3=1,3
      rmsg=rmsg+imove(i2)*d(i3,i2)*d(i3,i2)
      nmovh=nmovh+imove(i2)
      end do
      end do
      rmsg=sqrt(rmsg/float(nmovh*3))

      if (iopt.ne.1) write (81,*)rmsg
      caldat(ndata2)=rmsg
      end if

      if (isort(i1).eq.11) then
      if (iopt.ne.1) write (81,*)1000.0*pressu
      caldat(ndata2)=1000.0*pressu
      end if

   25 continue
      
   50 return
************************************************************************
*                                                                      *
*     Format part                                                      *
*                                                                      *
************************************************************************
  100 format (i3,a60)
  200 format (5i3)
      end
************************************************************************
************************************************************************

      subroutine readgeo

************************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkconst.blk"
#include "cbkdistan.blk"
#include "cbkinit.blk"
#include "cbkqa.blk"
#include "cbksrtbon1.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "small.blk"
      character*80 qromb
      character*25 qfileh
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readgeo'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

      if (ngeofor.eq.-1) return

********************************************************************** 
*                                                                    *
*     Read in system geometry                                        *
*                                                                    *
********************************************************************** 
      if (ngeofor.eq.0) then
      call readdelphi (qfileh,iend,naold)
      namov=na
      end if
 
      if (ngeofor.eq.1) then
      call readbgf(iend,naold)
      end if
 


      if (ngeofor.eq.2) then
**********************************************************************
*                                                                    *
*     Read in free format (xmol) geometry                            *
*                                                                    *
**********************************************************************
      qr='1'
      read (3,'(i6)')na
      namov=na
      read (3,'(a60)')qmol
      do i1=1,na
      read (3,'(a80)')qromb
      ifirstchar=80
      do i2=1,80
      if (qromb(i2:i2).ne.' '.and.i2.lt.ifirstchar) ifirstchar=i2
      end do
      read (qromb(ifirstchar:80),'(a2)')qa(i1)
      read (qromb(ifirstchar+2:80),*)c(i1,1),c(i1,2),c(i1,3)
      qlabel(i1)=qa(i1)
      qresi1(i1)='   '
      qresi2(i1)=' '
      qresi3(i1)='     '
      qffty(i1)='     '
      end do
      ibity=1
      axiss(1)=-1.0
      end if
 
 
      if (ngeofor.eq.3) then
**********************************************************************
*                                                                    *
*     Read in ChemDraw CC1-file                                      *
*                                                                    *
**********************************************************************
      qr='1'
      read (3,*)na
      namov=na
      read (3,'(a60)')qmol
      do i1=1,na
      read (3,'(2x,a2,5x,3f12.6)')qa(i1),c(i1,1),c(i1,2),c(i1,3)
      end do
      end if

      if (ngeofor.eq.4) then
**********************************************************************
*                                                                    *
*     Read in .pdb-format                                            *
*                                                                    *
**********************************************************************
      qr='C'
      call readpdb(iendf)
      namov=na
      ibity=1
      axiss(1)=-1.0
      qfile(nprob)=qmol
      if (iendf.eq.1) stop 'End-of-file while reading in .pdb'
      end if

**********************************************************************
*                                                                    *
*     Set up periodic system                                         *
*                                                                    *
**********************************************************************
      axis(1)=axiss(1)
      axis(2)=axiss(2)
      axis(3)=axiss(3)
      angle(1)=angles(1)
      angle(2)=angles(2)
      angle(3)=angles(3)
      if (axiss(1).lt.zero) then
      axis(1)=axis1
      axis(2)=axis2
      axis(3)=axis3
      angle(1)=90.0
      angle(2)=90.0
      angle(3)=90.0
      end if
      halfa=angle(1)*dgrrdn
      hbeta=angle(2)*dgrrdn
      hgamma=angle(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axis(1)*sinbet*sinphi
      tm21=axis(1)*sinbet*cosphi
      tm31=axis(1)*cosbet
      tm22=axis(2)*sinalf
      tm32=axis(2)*cosalf
      tm33=axis(3)

      return
      end
************************************************************************
************************************************************************

      subroutine readdelphi (qfileh,iend,naold)

************************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkconst.blk"
#include "cbkd.blk"
#include "cbkdcell.blk"
#include "cbkdistan.blk"
#include "cbkff.blk"
#include "cbkh.blk"
#include "cbkinit.blk"
#include "cbkqa.blk"
#include "cbkrestr.blk"
#include "cbksrtbon1.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
      character*25 qfileh
********************************************************************** 
*                                                                    *
*     Read in geometries in Delphi-format (xyz)                      *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readdelphi'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

      if (imodfile.eq.1) then
      open (3,file=qfileh,status='old')
      end if
      nmmax=nmmaxold
      nfc=nfcold
      ibity=1
      iredo=1
      endpo=endpoold
      icell=icellold
      icelo2=icelo2old
      iend=0
      read (3,1000,end=900)qr,qmol 
********************************************************************** 
*                                                                    *
*     Read in restraint information (optional)                       *
*                                                                    *
********************************************************************** 
      if (qr.eq.'R'.or.qr.eq.'P'.or.qr.eq.'X') then
      qmol=qmol(7:60)
      qmolset(nuge)=qmol
      read (18,1070,end=4,err=4) nrestra
      do i1=1,nrestra
      read (18,1090)irstra(i1,1),irstra(i1,2),rrstra(i1),vkrstr(i1),
     $vkrst2(i1),rrcha(i1)
      end do
    4 continue
      end if
********************************************************************** 
*                                                                    *
*     Read in torsion restraint information (optional)               *
*                                                                    *
********************************************************************** 
      if (qr.eq.'T'.or.qr.eq.'X') then
      if (qr.eq.'T') then
      qmol=qmol(7:60)
      qmolset(nuge)=qmol
      end if
      read (28,1070,end=6,err=6) nrestrat
      do i1=1,nrestrat
      read (28,1091)irstrat(i1,1),irstrat(i1,2),irstrat(i1,3),
     $irstrat(i1,4),trstra(i1),vkrt(i1),vkr2t(i1),rtcha(i1)
      end do
    6 continue
      end if
********************************************************************** 
*                                                                    *
*     Read in valency angle restraint information (optional)         *
*                                                                    *
********************************************************************** 
      if (qr.eq.'V') then
      qmol=qmol(7:60)
      qmolset(nuge)=qmol
      read (38,1070,end=7,err=7) nrestrav
      do i1=1,nrestrav
      read (38,1092)irstrav(i1,1),irstrav(i1,2),irstrav(i1,3),
     $vrstra(i1),vkrv(i1),vkr2v(i1)
      end do
    7 continue
      end if
********************************************************************** 
*                                                                    *
*     Read in geometry                                               *
*                                                                    *
********************************************************************** 
      ibh2=0
      iequi=1
      iexco=0
      if (nequi.gt.0) iequi=0
      axiss(1)=-1.0

      if (qr.eq.'O'.or.qr.eq.'L') stop 'Not xyz-format'

      if (qr.eq.'I') then      !Delphi internal coordinates
      if (nsurp.ge.2) stop 'Int.coordinates only with 1 gemetry' 
      call staint
      goto 20
      end if

      if (qr.eq.'B') then      !Previous geometry with volume reduction
      read (3,*)
      vred=(1.0-0.01*volcha)**(0.33333)
      iexco=1
      na=naold
      do i1=1,3
      qmol=qmol
      axiss(i1)=vred*axis(i1)
      angles(i1)=angle(i1)
      do i2=1,na
      c(i2,i1)=vred*c(i2,i1)
      end do
      end do

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      ibity=2

      goto 20
      end if

      if (qr.eq.'S') then      !Previous geometry with volume expansion
      read (3,*)
      vexp=(1.0+0.01*volcha)**(0.33333)
      na=naold
      iexco=1
      do i1=1,3
      qmol=qmol
      axiss(i1)=vexp*axis(i1)
      angles(i1)=angle(i1)
      do i2=1,na
      c(i2,i1)=vexp*c(i2,i1)
      end do
      end do

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      ibity=2

      goto 20
      end if

      if (qr.eq.'F'.or.qr.eq.'Y'.or.qr.eq.'3'.or.qr.eq.'5'.
     $or.qr.eq.'P') then 
      kx=0
      ky=0
      kz=0
      ibity=2
      read(3,1005)axiss(1),axiss(2),axiss(3)
      read(3,1005)angles(1),angles(2),angles(3)

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)

      end if

      if (qr.eq.'M'.or.qr.eq.'A') then 
      nmmsav=nmm
      nmm=2 
      end if

      if (qr.eq.'A') nmm=1

      if (qr.eq.'D') then
      endpo=endpo/25 
      nmmax=nmmax*5
      qruid='HIGH PRECISION'
      end if

      if (qr.eq.'H') then
      nmmax=nmmax/10
      qruid='LOW PRECISION'
      end if

      if (qr.eq.'1'.or.qr.eq.'5') then
      nmm=1
      nmmax=1
      qruid='SINGLE POINT'
      end if

      if (qr.eq.'Y') then
      icell=0
      qruid='NO CELL OPT'
      end if

   10 read (3,1100,end=20,err=20)ir,qa(na+1),(c(na+1,i2),i2=1,3)
      qlabel(na+1)=qa(na+1)
      qresi1(na+1)='   '
      qresi2(na+1)=' '
      qresi3(na+1)='     '
      qffty(na+1)='     '
      if (ir.eq.0) goto 20
      na=na+1

      if (na.gt.nat) then
      write (*,*)'Maximum number of atom exceeded ',na,nat
      stop 'Maximum number of atoms exceeded'
      end if

      goto 10
   20 continue

      if (imodfile.eq.1) close (3)

      return
  900 iend=1
      return
 1000 format (2x,a1,1x,a60)
 1005 format (3f10.4)
 1070 format (i3)
 1090 format (2i4,2f8.4,f8.6,f10.8)
 1091 format (4i4,2f8.4,3f8.6)
 1092 format (3i4,2f8.4,2f8.6)
 1100 format (i4,1x,a2,3x,3d22.15,1x,a5,1x,i5)
      end
************************************************************************
************************************************************************

      subroutine readbgf(iendf,naold)

************************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkcha.blk"
#include "cbkcharmol.blk"
#include "cbkconst.blk"
#include "cbkd.blk"
#include "cbkdcell.blk"
#include "cbkdistan.blk"
#include "cbkenergies.blk"
#include "cbkff.blk"
#include "cbkh.blk"
#include "cbkimove.blk"
#include "cbkinit.blk"
#include "cbkqa.blk"
#include "cbkrestr.blk"
#include "cbksrtbon1.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
      character*80 qromb
      character*2 qrom
      character*5 quen
      character*5 qlabhulp
      character*25 qfileh
      character*200 qhulp
********************************************************************** 
*                                                                    *
*     Read in BIOGRAF-geometry                                       *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readbgf'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

      iendf=0
      ienread=0
      iredo=0
      qremark(1)=' '
      enmol=zero
      formol=zero
c$$$      if (imodfile.eq.1) then
c$$$      open (3,file=qfileh,status='old')
c$$$      end if
      open (3,file='fort.3',status='old')
      read (3,'(a40)',end=900)qromb
      ibity=0
      if (qromb(1:6).eq.'BIOGRF') ibity=1
      if (qromb(1:6).eq.'XTLGRF') ibity=2
      if (ibity.eq.0) then
      write (*,*)qromb(1:6)
      stop 'Unknown Biograf-file'
      end if
      read (qromb,'(6x,i4)')ibgfversion
      if (ibity.eq.1) qr='C'
      if (ibity.eq.2) qr='F'
      iremark=0
      iformat=0
      iline=0
      iexco=0
      iruid=1
      vvol=1.0
      nmcharge=0
      nmmax=nmmaxold
      nfc=nfcold
      ncha=nchaold
      endpo=endpoold
      icell=icellold
      icelo2=icelo2old
      axiss(1)=-1.0

   30 read (3,'(a200)',end=46,err=40)qhulp
      qstrana1(1:200)=qhulp
      iline=iline+1
      irecog=0

      if (qhulp(1:6).eq.'DESCRP') then
      read (qhulp,'(7x,a40)',end=46,err=46)qmol
      irecog=1
      end if

      if (qhulp(1:6).eq.'REMARK') then
      if (iremark.lt.20) iremark=iremark+1
      read (qhulp,'(7x,a40)',end=46,err=46)qremark(iremark)
      irecog=1
      end if

      if (qhulp(1:6).eq.'FORMAT') then
      if (iformat.lt.20) iformat=iformat+1
      read(qhulp,'(7x,a40)',end=46,err=46)qformat(iformat)
      irecog=1
      end if

      if (qhulp(1:7).eq.'VCHANGE') then
      read (qhulp(8:60),*)vvol
      vred=(1.0+(vvol-1.0))**(0.33333333)
      iexco=1
      na=naold
      qmol=qmol
      do i1=1,3
      axiss(i1)=vred*axis(i1)
      angles(i1)=angle(i1)
      do i2=1,na
      cglobal(i2,i1)=vred*cglobal(i2,i1)
      end do
      end do

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      ibity=2
      irecog=1
      end if

      if (qhulp(1:7).eq.'VCHANGX') then
      read (qhulp(8:60),*)vvol
      vred=vvol
      iexco=1
      na=naold
      qmol=qmol
      do i1=1,3
      axiss(i1)=axis(i1)
      angles(i1)=angle(i1)
      do i2=1,na
      cglobal(i2,i1)=cglobal(i2,i1)
      end do
      end do

      axiss(1)=vred*axiss(1)
      do i2=1,na
      cglobal(i2,1)=vred*cglobal(i2,1)
      end do

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      ibity=2
      irecog=1
      end if

      if (qhulp(1:7).eq.'VCHANGY') then
      read (qhulp(8:60),*)vvol
      vred=vvol
      iexco=1
      na=naold
      qmol=qmol
      do i1=1,3
      axiss(i1)=axis(i1)
      angles(i1)=angle(i1)
      do i2=1,na
      cglobal(i2,i1)=cglobal(i2,i1)
      end do
      end do

      axiss(2)=vred*axiss(2)
      do i2=1,na
      cglobal(i2,2)=vred*cglobal(i2,2)
      end do

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      ibity=2
      irecog=1
      end if

      if (qhulp(1:7).eq.'VCHANGZ') then
      read (qhulp(8:60),*)vvol
      vred=vvol
      iexco=1
      na=naold
      qmol=qmol

      do i1=1,3
      axiss(i1)=axis(i1)
      angles(i1)=angle(i1)
      do i2=1,na
      cglobal(i2,i1)=cglobal(i2,i1)
      end do
      end do

      axiss(3)=vred*axiss(3)
      do i2=1,na
      cglobal(i2,3)=vred*cglobal(i2,3)
      end do

      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      ibity=2
      irecog=1
      end if

      if (qhulp(1:6).eq.'CRYSTX') then
      read (qhulp,'(8x,6f11.5)',end=46,err=46)axiss(1),
     $axiss(2),axiss(3),angles(1),angles(2),angles(3)
      kx=0
      ky=0
      kz=0
      halfa=angles(1)*dgrrdn
      hbeta=angles(2)*dgrrdn
      hgamma=angles(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axiss(1)*sinbet*sinphi
      tm21=axiss(1)*sinbet*cosphi
      tm31=axiss(1)*cosbet
      tm22=axiss(2)*sinalf
      tm32=axiss(2)*cosalf
      tm33=axiss(3)
      kx=int(2.0*swb/tm11)
      ky=int(2.0*swb/tm22)
      kz=int(2.0*swb/tm33)
      qr='F'
      if (nmmax.eq.1.and.nmmaxold.gt.1) qr='5'
      if (icell.eq.0.and.icellold.gt.0) qr='Y'
      ibity=2
      irecog=1
      end if

      if (qhulp(1:6).eq.'PERIOD') then
      read (qhulp,'(7x,i3)',end=46,err=46)iperiod
      irecog=1
      end if

      if (qhulp(1:4).eq.'AXES') then
      read (qhulp,'(7x,a3)',end=46,err=46)qbgfaxes
      irecog=1
      end if

      if (qhulp(1:6).eq.'SGNAME') then
      read (qhulp,'(7x,a3)',end=46,err=46)qbgfsgn
      irecog=1
      end if

*     if (qhulp(1:5).eq.'CELLS') then
*     read (qhulp,'(7x,*)',end=40,err=40)kx,ky,kz
*     irecog=1
*     end if

      if (qhulp(1:6).eq.'HETATM') then
      if (ibgfversion.lt.400) then
      read (qhulp,
     $'(7x,i5,1x,a5,1x,a3,1x,a1,1x,a5,3f10.5,1x,a5,i3,i2,1x,f8.5)'
     $,end=40,err=40)
     $ir,qlabel(na+1),qresi1(na+1),qresi2(na+1),qresi3(na+1),
     $cglobal(na+1,1),cglobal(na+1,2),
     $cglobal(na+1,3),qffty(na+1),ibgr1(na+1),ibgr2(na+1),
     $chgglobal(na+1)
      else
      stop 'Unsupported Biograf-version'
      end if
      qlabhulp=qlabel(na+1)
      if (qlabhulp(1:1).eq.' ') qlabhulp=qlabhulp(2:5)
      if (qlabhulp(1:1).eq.' ') qlabhulp=qlabhulp(2:4)
      if (qlabhulp(1:1).eq.' ') qlabhulp=qlabhulp(2:3)
      if (qlabhulp(1:1).eq.'C ') qa(na+1)='C '
      if (qlabhulp(1:2).eq.'Ca') qa(na+1)='Ca'
      if (qlabhulp(1:2).eq.'Cl') qa(na+1)='Cl'
      if (qlabhulp(1:2).eq.'Cu') qa(na+1)='Cu'
      if (qlabhulp(1:2).eq.'Co') qa(na+1)='Co'
      if (qlabhulp(1:1).eq.'H ') qa(na+1)='H '
      if (qlabhulp(1:2).eq.'He') qa(na+1)='He'
      if (qlabhulp(1:1).eq.'N ') qa(na+1)='N '
      if (qlabhulp(1:2).eq.'Ni') qa(na+1)='Ni'
      if (qlabhulp(1:1).eq.'O ') qa(na+1)='O '
      if (qlabhulp(1:1).eq.'B ') qa(na+1)='B '
      if (qlabhulp(1:1).eq.'F ') qa(na+1)='F '
      if (qlabhulp(1:2).eq.'Fe') qa(na+1)='Fe'
      if (qlabhulp(1:1).eq.'P ') qa(na+1)='P '
      if (qlabhulp(1:1).eq.'S ') qa(na+1)='S '
      if (qlabhulp(1:1).eq.'Y ') qa(na+1)='Y '
      if (qlabhulp(1:2).eq.'Al ') qa(na+1)='Al'
      if (qlabhulp(1:2).eq.'Au ') qa(na+1)='Au'
      if (qlabhulp(1:2).eq.'Si') qa(na+1)='Si'
      if (qlabhulp(1:2).eq.'Pt') qa(na+1)='Pt'
      if (qlabhulp(1:2).eq.'Mo') qa(na+1)='Mo'
      if (qlabhulp(1:2).eq.'Mg') qa(na+1)='Mg'
      if (qlabhulp(1:2).eq.'Ar') qa(na+1)='Ar'
      if (qlabhulp(1:2).eq.'Zr') qa(na+1)='Zr'
      if (qlabhulp(1:2).eq.'Ti') qa(na+1)='Ti'
      if (qlabhulp(1:2).eq.'Ru') qa(na+1)='Ru'
      if (qlabhulp(1:2).eq.'Ba') qa(na+1)='Ba'
      if (qlabhulp(1:2).eq.'Bi') qa(na+1)='Bi'
      if (qlabhulp(1:2).eq.'Li') qa(na+1)='Li'
      if (qlabhulp(1:2).eq.'V ') qa(na+1)='V '
      if (qlabhulp(1:2).eq.'X ') qa(na+1)='X '
      na=na+1
      if (na.gt.nattot) then
      write (*,*)'Number of atoms:read ',na
      write (*,*)'Maximum number of atoms: ',nattot
      stop 
     $'Maximum number of atoms exceeded; increase nattot in cbka.blk'
      end if
      irecog=1
      end if

      if (qhulp(1:6).eq.'RUTYPE') then          !run-type identifiers
      irecrun=0
      read (qhulp,'(7x,a40)',end=46,err=46)qruid

      if (qruid(1:10).eq.'NORMAL RUN') then
      iruid=0
      irecrun=1
      end if

      if (qruid(1:14).eq.'HIGH PRECISION') then
      endpo=endpo/25 
      nmmax=nmmax*5
      qr='D'
      iruid=1
      irecrun=1
      end if

      if (qruid(1:13).eq.'LOW PRECISION') then
      nmmax=nmmax/10
      qr='H'
      iruid=1
      irecrun=1
      end if

      if (qruid(1:12).eq.'SINGLE POINT') then
      iruid=1
      nmmax=1
      qr='1'
      if (ibity.eq.2) qr='5'
      irecrun=1
      end if

      if (qruid(1:11).eq.'NO CELL OPT') then
      iruid=1
      icell=0
      if (ibity.eq.2) qr='Y'
      irecrun=1
      end if

      if (qruid(1:8).eq.'CELL OPT') then
      iruid=1
      icell=1
      iexco=0   !Override from VCHANGE
      read (qruid,'(8x,i6)',end=46,err=46)ncellopt
      if (ncellopt.eq.2) icell=2 !cell optimisation during energy minimisation 
      if (ncellopt.eq.3) icelo2=4  !c/a optimisation
      if (ncellopt.eq.4) icelo2=1  !only a optimisation
      if (ncellopt.eq.5) icelo2=2  !only b optimisation
      if (ncellopt.eq.6) icelo2=3  !only c optimisation
      if (ncellopt.eq.7) then
      icelo2=4  !c/a optimisation
      icell=2 !cell optimisation during energy minimisation
      end if
      if (ibity.eq.2) qr='F'
      irecrun=1
      end if

      if (qruid(1:6).eq.'MAXMOV') then
      iruid=1
      read (qruid,'(6x,i6)',end=46,err=46)nfc
      irecrun=1
      end if

      if (qruid(1:4).eq.'REDO') then
      iruid=1
      read (qruid,'(4x,i6)',end=46,err=46)iredo
      irecrun=1
      end if

      if (qruid(1:5).eq.'MAXIT') then
      iruid=1
      read (qruid,'(6x,i6)',end=46,err=46)nmmax
      if (qruid(14:18).eq.'ENDPO') then
      read (qruid,'(18x,f6.3)',end=46,err=46)endpo
      end if
      irecrun=1
      end if
      if (qruid(1:5).eq.'ENDPO') then
      iruid=1
      read (qruid,'(6x,f6.3)',end=46,err=46)endpo
      irecrun=1
      end if
      
      if (qruid(1:9).eq.'CHARGEMET') then
      iruid=1
      read (qruid,'(9x,i6)',end=46,err=46)ncha
      irecrun=1
      end if

      if (irecrun.eq.0) then
      write (*,*)'Warning: ignored RUTYPE identifier ',qruid(1:12)
      end if

      irecog=1
      end if

      if (qhulp(1:14).eq.'BOND RESTRAINT') then  
      nrestra=nrestra+1
      istart=15
      call stranal(istart,iend,vout,iout,1)
      irstra(nrestra,1)=iout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      irstra(nrestra,2)=iout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      rrstra(nrestra)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      vkrstr(nrestra)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      vkrst2(nrestra)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      rrcha(nrestra)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      itstart(nrestra)=iout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      itend(nrestra)=iout
      istart=iend
*     read (qhulp,'(15x,2i4,f8.4,f8.2,f8.5,f10.7)',end=46,err=46)
*    $irstra(nrestra,1),irstra(nrestra,2),rrstra(nrestra),
*    $vkrstr(nrestra),vkrst2(nrestra),rrcha(nrestra)
      qr='R'
      irecog=1
      end if

      if (qhulp(1:15).eq.'ANGLE RESTRAINT') then  
      nrestrav=nrestrav+1
      read (qhulp,'(16x,3i4,2f8.2,f8.4,f9.6)',end=46,err=46)
     $irstrav(nrestrav,1),irstrav(nrestrav,2),irstrav(nrestrav,3),
     $vrstra(nrestrav),vkrv(nrestrav),vkr2v(nrestrav),
     $rvcha(nrestrav)
      qr='V'
      irecog=1
      end if

      if (qhulp(1:17).eq.'TORSION RESTRAINT') then  
      nrestrat=nrestrat+1
      read (qhulp,'(18x,4i4,2f8.2,f8.4,f9.6)',end=46,err=46)
     $irstrat(nrestrat,1),irstrat(nrestrat,2),irstrat(nrestrat,3),
     $irstrat(nrestrat,4),trstra(nrestrat),vkrt(nrestrat),
     $vkr2t(nrestrat),rtcha(nrestrat)
      qr='T'
      irecog=1
      end if

      if (qhulp(1:16).eq.'MASCEN RESTRAINT') then  
      nrestram=nrestram+1
      istart=17
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      irstram(nrestram,1)=0
      if (qstrana2.eq.'x') irstram(nrestram,1)=1 
      if (qstrana2.eq.'y') irstram(nrestram,1)=2 
      if (qstrana2.eq.'z') irstram(nrestram,1)=3 
      if (qstrana2.eq.'p') irstram(nrestram,1)=4 !fixed center of mass
      if (irstram(nrestram,1).eq.0) 
     $stop 'Error in mass centre restraint'
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      irstram(nrestram,2)=iout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      irstram(nrestram,3)=iout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      rmstra1(nrestram)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      if (irstram(nrestram,1).le.3) irstram(nrestram,4)=iout
      if (irstram(nrestram,1).eq.4) rmstra2(nrestram)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      if (irstram(nrestram,1).le.3) irstram(nrestram,5)=iout
      if (irstram(nrestram,1).eq.4) rmstra3(nrestram)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      if (irstram(nrestram,1).le.3) rmstra2(nrestram)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      if (irstram(nrestram,1).le.3) rmstra3(nrestram)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      if (irstram(nrestram,1).le.3) rmcha(nrestram)=vout
      irecog=1
      end if

      if (qhulp(1:9).eq.'MOLCHARGE') then
      nmcharge=nmcharge+1
      istart=10
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      iat1mc(nmcharge)=iout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      iat2mc(nmcharge)=iout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      vmcha(nmcharge)=vout
      irecog=1
      end if
      
      if (qhulp(1:8).eq.'FIXATOMS') then
      istart=9
      call stranal(istart,iend,vout,iout,1)
      if1=iout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      if2=iout
      do i12=if1,if2
      imove(i12)=0
      end do
      irecog=1
      end if

      if (qhulp(1:11).eq.'UNIT ENERGY') then
      eenconv=zero
      read (qhulp,'(14x,a5)',end=46,err=46)quen
      if (quen.eq.'eV') eenconv=23.0408
      if (quen.eq.'EV') eenconv=23.0408
      if (quen.eq.'ev') eenconv=23.0408
      if (quen.eq.'h') eenconv=627.5
      if (quen.eq.'H') eenconv=627.5
      if (quen.eq.'kcal') eenconv=1.0
      if (quen.eq.'kCal') eenconv=1.0
      if (quen.eq.'KCAL') eenconv=1.0
      if (eenconv.eq.zero) then 
      write (*,*)quen,': unknown energy unit; assuming kcal/mol'
      eenconv=1.0
      end if
      irecog=1
      end if

      if (qhulp(1:6).eq.'ENERGY') then
      read (qhulp(7:80),*,end=46,err=46)enmol
      ienread=1
      irecog=1
      end if

      if (qhulp(1:6).eq.'GEOUPD') then
      icgeopt(nprob)=0
      icgeo=0
      irecog=1
      end if
 
      if (qhulp(1:9).eq.'NO GEOUPD') then
      icgeopt(nprob)=1
      icgeo=1
      irecog=1
      end if

      if (qhulp(1:9).eq.'FREQUENCY') then
      ifreqset(nprob)=1
      ifreq=1
      irecog=1
      end if

*     if (qhulp(1:5).eq.'FORCE') then
*     read (qhulp(6:80),*,end=46,err=46)formol
*     ienread=1
*     irecog=1
*     end if

      if (qhulp(1:6).eq.'FFIELD') goto 30 
      if (qhulp(1:6).eq.'CONECT') goto 30
      if (qhulp(1:5).eq.'ORDER') goto 30
      if (qhulp(1:1).eq.'#') goto 30
      if (qhulp(1:3).eq.'END') goto 45

      if (irecog.eq.0) then
      write (*,*)'Warning: ignored line starting with: ',qhulp(1:10)
      end if
 
      goto 30
      
   40 write (*,*)'Error on line ',iline+1,' of Biograf-input'
      stop
   45 read (3,*,err=46,end=46)
   46 continue
      if (ienread.eq.1) then 
      if (eenconv.eq.zero) then
      write (*,*)'No energy unit given; assuming kcal/mol' 
      eenconv=1.0
      end if
      enmol=enmol*eenconv                !Convert energies to kcal/mol
      end if
      
      namov=0                            !calculate number of moving atoms
      do i1=1,na
      if (imove(i1).eq.1) namov=namov+1
      end do

      if (imodfile.eq.1) close (3)
      return
  900 iendf=1
      return
      end
************************************************************************
************************************************************************

      subroutine readpdb (iendf)

************************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkqa.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
#include "cbksrtbon1.blk"
      character*200 qhulp
********************************************************************** 
*                                                                    *
*     Read in .pdb-geometry                                          *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readpdb'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      iendf=1
      qmol='pdb_in'
    5 read (3,'(a200)',end=10,err=900) qhulp
      qstrana1(1:200)=qhulp
      istart=1
      call stranal(istart,iend,vout,iout,1)
      istart=iend

      if (qstrana2(1:6).eq.'HEADER') then
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      qmol=qstrana2(1:20)
      end if

      if (qstrana2(1:6).eq.'HETATM'.or.qstrana2(1:4).eq.'ATOM') then
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      qa(na+1)=qstrana2(1:2)
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      c(na+1,1)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      c(na+1,2)=vout
      call stranal(istart,iend,vout,iout,1)
      istart=iend
      c(na+1,3)=vout
      na=na+1
      end if

      if (qstrana2(1:3).eq.'END'.or.qstrana2(2:4).eq.'END') then
      iendf=0
      goto 10
      end if

      goto 5
   10 continue
      return
  900 write (*,*)'Error reading in .pdb-format'
      stop 'Error reading in .pdb-format'
      end 
************************************************************************
************************************************************************

      subroutine readtreg

************************************************************************
#include "cbka.blk"
#include "cbktregime.blk"
#include "control.blk"
#include "small.blk"
      dimension isumattreg(mtreg)
      character*200 qrom
********************************************************************** 
*                                                                    *
*     Read in temperature regime                                     *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readtreg'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      ntrc=0
      open (19,file='tregime.in',status='old',err=60)
   10 read (19,'(a200)',end=50,err=900)qrom
      qstrana1(1:200)=qrom
      if (qrom(1:1).eq.'#') goto 10
      istart=1
      ntrc=ntrc+1
      if (ntrc.gt.mtreg) then
      write (*,*)'Too many temperature regimes in tregime.in;',
     $' inrease mtreg in cbka.blk'
      stop 'Too many temperature regimes in tregime.in'
      end if
      call stranal(istart,iend,vout,iout,1)
      nittc(ntrc)=iout
      istart=iend

      if (ntrc.gt.1) then
      if (nittc(ntrc).lt.nittc(ntrc-1)) then
      ntrc=ntrc-1
      write (*,*)'Warning: wrong order or empty line in tregime.in'
      write (*,*)'Ignored lines below iteration:',nittc(ntrc)
      goto 50
      end if
      end if

      call stranal(istart,iend,vout,iout,1)
      nntreg(ntrc)=iout
      if (nntreg(ntrc).gt.mtzone) then
      write (*,*)'Too many temperature zones in tregime.in;',
     $' inrease mtzone in cbka.blk'
      stop 'Too many temperature zones in tregime.in'
      end if
      istart=iend
      isumattreg(ntrc)=0
      do i1=1,nntreg(ntrc)
      call stranal(istart,iend,vout,iout,1)
      ia1treg(ntrc,i1)=iout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      ia2treg(ntrc,i1)=iout
      istart=iend
      isumattreg(ntrc)=isumattreg(ntrc)+1+ia2treg(ntrc,i1)-
     $ia1treg(ntrc,i1)
      call stranal(istart,iend,vout,iout,1)
      tsettreg(ntrc,i1)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      tdamptreg(ntrc,i1)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      dttreg(ntrc,i1)=vout
      istart=iend
      end do
      goto 10
   50 continue
      close (19)
   60 continue
**********************************************************************
*                                                                    *
*     Check consistency temperature programs in tregime.in           *
*                                                                    *
**********************************************************************
      if (ntrc.gt.0) then
      do i1=1,ntrc
      if (isumattreg(i1).ne.na) then
      write (*,*)'Inconsistency in temperature regime nr.',i1
      write (*,*)'Number of atoms defined in tregime.in:',
     $isumattreg(i1)
      write (*,*)'Number of atoms in system:',na
      stop 'Inconsistency in tregime.in'
      end if
      end do
      end if

      return
  900 stop 'Error reading tregime.in'
      end
************************************************************************
************************************************************************

      subroutine readvreg

************************************************************************
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkvregime.blk"
#include "control.blk"
      character*200 qrom
********************************************************************** 
*                                                                    *
*     Read in volume regime                                          *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readvreg'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      nvrc=0
      open (19,file='vregime.in',status='old',err=60)
   10 read (19,'(a200)',end=50,err=900)qrom
      qstrana1(1:200)=qrom
      if (qrom(1:1).eq.'#') goto 10
      istart=1
      nvrc=nvrc+1
      if (nvrc.gt.mvreg) then
      write (*,*)'Too many volume regimes in vregime.in;',
     $' inrease mvreg in cbka.blk'
      stop 'Too many volume regimes in vregime.in'
      end if

      call stranal(istart,iend,vout,iout,1)
      nitvc(nvrc)=iout
      istart=iend

      if (nvrc.gt.1) then
      if (nitvc(nvrc).lt.nitvc(nvrc-1)) then
      nvrc=nvrc-1
      write (*,*)'Warning: wrong order or empty line in vregime.in'
      write (*,*)'Ignored lines below iteration:',nitvc(nvrc)
      goto 50
      end if
      end if

      call stranal(istart,iend,vout,iout,1)
      nnvreg(nvrc)=iout
      if (nnvreg(nvrc).gt.mvzone) then
      write (*,*)'Too many volume regimes in vregime.in;',
     $' inrease mvzone in cbka.blk'
      stop 'Too many volume zones in vregime.in'
      end if
      istart=iend
      do i1=1,nnvreg(nvrc)
      call stranal(istart,iend,vout,iout,1)
      if (qstrana2(1:1).ne.'a'.and.qstrana2(1:1).ne.'b'.and.
     $qstrana2(1:1).ne.'c'.and.qstrana2(1:4).ne.'alfa'.and.
     $qstrana2(1:4).ne.'beta'.and.qstrana2(1:5).ne.'gamma') then
      write (*,*)qstrana2
      write (*,*)'Invalid cell parameter type in vregime.in ;',
     $' use a,b,c,alfa,beta or gamma'
      stop 'Invalid cell parameter type in vregime.in'
      end if
      qvtype(nvrc,i1)=qstrana2
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      dvvreg(nvrc,i1)=vout
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      ivsca(nvrc,i1)=1
      if (qstrana2(1:1).eq.'n') ivsca(nvrc,i1)=0
      istart=iend
      end do
      goto 10
   50 continue
      close (19)
   60 continue
      return
  900 stop 'Error reading vregime.in'
      end
************************************************************************
************************************************************************

      subroutine readereg

************************************************************************
#include "cbka.blk"
#include "cbkeregime.blk"
#include "control.blk"
      character*200 qrom
********************************************************************** 
*                                                                    *
*     Read in electric field regime                                  *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readereg'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      nerc=0
      open (19,file='eregime.in',status='old',err=60)
   10 read (19,'(a200)',end=50,err=900)qrom
      qstrana1(1:200)=qrom
      if (qrom(1:1).eq.'#') goto 10
      istart=1
      nerc=nerc+1
      if (nerc.gt.mereg) then
      write (*,*)'Too many electric field regimes in eregime.in;',
     $' inrease mereg in cbka.blk'
      stop 'Too many electric field regimes in eregime.in'
      end if
      call stranal(istart,iend,vout,iout,1)
      nitec(nerc)=iout

      if (nerc.gt.1) then
      if (nitec(nerc).lt.nitec(nerc-1)) then
      nerc=nerc-1
      write (*,*)'Warning: wrong order or empty line in eregime.in'
      write (*,*)'Ignored lines below iteration:',nitec(nerc)
      goto 50
      end if
      end if
      
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      nnereg(nerc)=iout
      if (nnereg(nerc).gt.mezone) then
      write (*,*)'Too many electric field zones in eregime.in;',
     $' inrease mezone in cbka.blk'
      stop 'Too many electric field zones in vregime.in'
      end if
      istart=iend
      do i1=1,nnereg(nerc)
      call stranal(istart,iend,vout,iout,1)
      if (qstrana2(1:1).ne.'x'.and.qstrana2(1:1).ne.'y'.and.
     $qstrana2(1:1).ne.'z') then
      write (*,*)qstrana2
      write (*,*)'Invalid field direction in eregime.in ;',
     $' use x,y or z'
      stop 'Invalid field direction in eregime.in'
      end if
      qetype(nerc,i1)=qstrana2
      istart=iend
      call stranal(istart,iend,vout,iout,1)
      ereg(nerc,i1)=vout
      istart=iend
      end do
      goto 10
   50 continue
      close (19)
   60 continue
      return
  900 stop 'Error reading vregime.in'
      end
************************************************************************
************************************************************************

      subroutine readaddmol

************************************************************************
#include "cbka.blk"
#include "cbkatomcoord.blk"
#include "cbkc.blk"
#include "cbkff.blk"
#include "cbkh.blk"
#include "control.blk"
      character*80 qromb
      character*200 qhulp
      character*5 qlabhulp
********************************************************************** 
*                                                                    *
*     Read in molecule coordinates. This molecule will be added to   *
*     the system at regular intervals                                *
*     Accepts only .bgf-format                                       *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readaddmol'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
********************************************************************** 
*                                                                    *
*     Set default values                                             *
*                                                                    *
********************************************************************** 
      iaddfreq=-1   !frequency of molecule addition; <0: no addition
      xadd=-9000.0  !x-coordinate for added molecule; <-5000.0: random
      yadd=-9000.0  !y-coordinate for added molecule; <-5000.0: random
      zadd=-9000.0  !z-coordinate for added molecule; <-5000.0: random
      iveladd=1     !1: random initial velocities; 2: read in velocities
                    !from addmol.vel
      addist=-1.00  !Minimum distance between added molecule and rest
                    !of system. < 0.0: do not check
      nadattempt=10  !Number of attempts at adding the molecule
      taddmol=-1.0  !Temperature added molecule. <0.0: system temperature
      open (19,file='addmol.bgf',status='old',err=60)
      read (19,'(a40)',end=900,err=900)qromb
      if (qromb(1:6).ne.'BIOGRF') then
      write (*,*)'addmol.bgf should start with BIOGRF'
      stop 'addmol.bgf should start with BIOGRF'
      end if
      naa=0
      iline=0
   30 read (19,'(a200)',end=900,err=900)qhulp
      irecog=0
      iline=iline+1

      if (qhulp(1:6).eq.'DESCRP') then
      irecog=1
      end if
 
      if (qhulp(1:6).eq.'FORMAT') then
      irecog=1
      end if

      if (qhulp(1:6).eq.'REMARK') then
      irecog=1
      end if

      if (qhulp(1:6).eq.'HETATM') then
      irecog=1
      read (qhulp,'(7x,i5,1x,a5,1x,3x,1x,1x,1x,5x,3f10.5)'
     $,end=900,err=900)
     $ir,qlabhulp,cadd(naa+1,1),cadd(naa+1,2),cadd(naa+1,3)
      if (qlabhulp(1:1).eq.' ') qlabhulp=qlabhulp(2:5)
      if (qlabhulp(1:1).eq.' ') qlabhulp=qlabhulp(2:4)
      if (qlabhulp(1:1).eq.' ') qlabhulp=qlabhulp(2:3)
      if (qlabhulp(1:1).eq.'C ') qadd(naa+1)='C '
      if (qlabhulp(1:2).eq.'Ca') qadd(naa+1)='Ca'
      if (qlabhulp(1:2).eq.'Cl') qadd(naa+1)='Cl'
      if (qlabhulp(1:2).eq.'Cu') qadd(naa+1)='Cu'
      if (qlabhulp(1:2).eq.'Co') qadd(naa+1)='Co'
      if (qlabhulp(1:1).eq.'H ') qadd(naa+1)='H '
      if (qlabhulp(1:2).eq.'He') qadd(naa+1)='He'
      if (qlabhulp(1:1).eq.'N ') qadd(naa+1)='N '
      if (qlabhulp(1:2).eq.'Ni') qadd(naa+1)='Ni'
      if (qlabhulp(1:1).eq.'O ') qadd(naa+1)='O '
      if (qlabhulp(1:1).eq.'B ') qadd(naa+1)='B '
      if (qlabhulp(1:1).eq.'F ') qadd(naa+1)='F '
      if (qlabhulp(1:2).eq.'Fe') qadd(naa+1)='Fe'
      if (qlabhulp(1:1).eq.'P ') qadd(naa+1)='P '
      if (qlabhulp(1:1).eq.'S ') qadd(naa+1)='S '
      if (qlabhulp(1:1).eq.'Y ') qadd(naa+1)='Y '
      if (qlabhulp(1:2).eq.'Al') qadd(naa+1)='Al'
      if (qlabhulp(1:2).eq.'Au') qadd(naa+1)='Au'
      if (qlabhulp(1:2).eq.'Si') qadd(naa+1)='Si'
      if (qlabhulp(1:2).eq.'Pt') qadd(naa+1)='Pt'
      if (qlabhulp(1:2).eq.'Mo') qadd(naa+1)='Mo'
      if (qlabhulp(1:2).eq.'Mg') qadd(naa+1)='Mg'
      if (qlabhulp(1:2).eq.'Ar') qadd(naa+1)='Ar'
      if (qlabhulp(1:2).eq.'Zr') qadd(naa+1)='Zr'
      if (qlabhulp(1:2).eq.'Ba') qadd(naa+1)='Ba'
      if (qlabhulp(1:2).eq.'X ') qadd(naa+1)='X '
      ityadd(naa+1)=0
      do i1=1,nso  !Find force field type
      if (qadd(naa+1).eq.qas(i1)) ityadd(naa+1)=i1 
      end do
      if (ityadd(naa+1).eq.0) then
      write (*,*) 'Unknown atom type:',qadd(naa+1)
      stop 'Unknown atom type'
      end if
      naa=naa+1
      end if

      if (qhulp(1:7).eq.'FREQADD') then
      irecog=1
      read (qhulp,'(8x,i6)',end=900,err=900) iaddfreq
      end if

      if (qhulp(1:6).eq.'VELADD') then
      irecog=1
      read (qhulp,'(8x,i6)',end=900,err=900) iveladd
      end if

      if (qhulp(1:6).eq.'STARTX') then
      irecog=1
      read (qhulp,'(7x,f8.2)',end=900,err=900) xadd
      end if

      if (qhulp(1:6).eq.'STARTY') then
      irecog=1
      read (qhulp,'(7x,f8.2)',end=900,err=900) yadd
      end if

      if (qhulp(1:6).eq.'STARTZ') then
      irecog=1
      read (qhulp,'(7x,f8.2)',end=900,err=900) zadd
      end if

      if (qhulp(1:6).eq.'ADDIST') then
      irecog=1
      read (qhulp,'(7x,f8.2)',end=900,err=900) addist
      end if

      if (qhulp(1:8).eq.'NATTEMPT') then
      irecog=1
      read (qhulp,'(9x,i6)',end=900,err=900) nadattempt
      end if

      if (qhulp(1:7).eq.'TADDMOL') then
      irecog=1
      read (qhulp,'(8x,f8.2)',end=900,err=900) taddmol
      end if

      if (qhulp(1:6).eq.'FFIELD') goto 30
      if (qhulp(1:6).eq.'CONECT') goto 30
      if (qhulp(1:5).eq.'ORDER') goto 30
      if (qhulp(1:1).eq.'#') goto 30
      if (qhulp(1:3).eq.'END') goto 45

      if (irecog.eq.0) then
      write (*,*)'Warning: ignored line starting with: ',qhulp(1:10)
      end if

      goto 30

   45 continue
      close (19)
      if (iveladd.eq.2) then
      open (19,file='addmol.vel',status='old',err=800)
      read (19,*)
      read (19,'(3d24.15)',err=850,end=850)
     $((veladd(j,i),j=1,3),i=1,naa)       
      close (19)
      end if
************************************************************************
*                                                                      *
*     Place molecule at origin                                         *
*                                                                      *
************************************************************************
      ccx=0.0
      ccy=0.0
      ccz=0.0
      do i1=1,naa
      ccx=ccx+cadd(i1,1)/float(naa)
      ccy=ccy+cadd(i1,2)/float(naa)
      ccz=ccz+cadd(i1,3)/float(naa)
      end do
      do i1=1,naa
      cadd(i1,1)=cadd(i1,1)-ccx
      cadd(i1,2)=cadd(i1,2)-ccy
      cadd(i1,3)=cadd(i1,3)-ccz
      end do

   60 continue
      return
  800 stop 'Error opening addmol.vel'
  850 stop 'Error or end of file reading addmol.vel'
  900 write (*,*)'Error or end-of-file reading addmol.bgf on line:',
     $iline
      return
      end
************************************************************************
**********************************************************************

      subroutine writegeo(nunit1)

**********************************************************************   
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkconst.blk"
#include "cbkqa.blk"
#include "cbkrestr.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
#include "cbksrtbon1.blk"
#include "cbkinit.blk"
**********************************************************************   
*                                                                    *
*     Copy new geometries to unit nunit1                             *
*                                                                    *
**********************************************************************   
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In writegeo'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      if (axiss(1).lt.zero) then
      if (nrestra.eq.0.and.nrestrat.eq.0.and.
     $nrestrav.eq.0) 
     $write (nunit1,300)qr,qmol
      if (nrestra.gt.0) write (nunit1,301)qr,
     $rrstra(1),qmol
      if (nrestrav.gt.0) write (nunit1,301)qr,
     $vrstra(1),qmol
      if (nrestrat.gt.0) write (nunit1,301)qr,
     $trstra(1),qmol
      else
      write (nunit1,310)qr,qmol
      write (nunit1,320)axiss(1),axiss(2),axiss(3)
      write (nunit1,320)angles(1),angles(2),angles(3)
      end if
      do i1=1,na
      if (nbiolab.ne.1) write (nunit1,400)i1,qa(i1),(c(i1,i2),i2=1,3)
      if (nbiolab.eq.1) write (nunit1,401)i1,qa(i1),(c(i1,i2),i2=1,3) !Delphi-format
      end do
      if (nbiolab.ne.1) write (nunit1,*)
     
      return

  300 format (2x,a1,1x,a60)
  301 format (2x,a1,1x,f6.2,a60)
  310 format (2x,a1,1x,a60)
  320 format (3f10.4)
  400 format (i4,1x,a2,3x,3(d21.14,1x),1x,a5,1x,i5)
  401 format (i3,2x,a2,3x,3(d21.14,1x),1x,a5,1x,i5)
      end
**********************************************************************
**********************************************************************

      subroutine writebgf(nunit1)

**********************************************************************   
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkcha.blk"
#include "cbkcharmol.blk"
#include "cbkconst.blk"
#include "cbkenergies.blk"
#include "cbkia.blk"
#include "cbkimove.blk"
#include "cbkinit.blk"
#include "cbkqa.blk"
#include "cbkrestr.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "opt.blk"
#include "cbksrtbon1.blk"
#include "small.blk"

      dimension qdir(3) 
      character*2 qt
      character*1 qdir
**********************************************************************   
*                                                                    *
*     Copy new Biograf-geometries to unit nunit1                     *
*                                                                    *
**********************************************************************   
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In newbgf'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
      irom=1
      qdir(1)='x'
      qdir(2)='y'
      qdir(3)='z'
      ibgfversion=200
      if (ibity.eq.1) write (nunit1,1500)ibgfversion
      if (ibity.eq.2) write (nunit1,1600)ibgfversion
*     if (qr.ne.'F'.and.qr.ne.'5'.and.qr.ne.'Y') 
*    $write (nunit1,1500)ibgfversion
*     if (qr.eq.'F'.or.qr.eq.'5'.or.qr.eq.'Y') 
*    $write (nunit1,1600)ibgfversion
      write (nunit1,1700)qmol
*     write (nunit1,1700)qkeyw(nprob)
      do i1=1,iremark
      write (nunit1,1800)qremark(i1)
      end do
      qruid='NORMAL RUN'
      if (iruid.eq.0) then
      write (nunit1,2000)
      else
      if (abs(endpo-endpoold).gt.1e-5) write (nunit1,2010)endpo
      if (nmmax.ne.nmmaxold) write (nunit1,2020)nmmax
      if (nfc.ne.nfcold) write (nunit1,2030)nfc
      if (ncha.ne.nchaold) write (nunit1,2036)ncha
      if (iredo.gt.1) write (nunit1,2035)iredo
      if (icell.ne.icellold) then
      if (icell.eq.0) write (nunit1,2033)
      if (icell.gt.0) write (nunit1,2034)ncellopt
      end if
      end if
      if (iexco.ne.0.and.nsurp.gt.0) then
      write (nunit1,2040)vvol
      write (nunit1,3500)
      write (nunit1,*)
      return
      end if
      if (nmcharge.gt.0) then
      do i3=1,nmcharge
      write (nunit1,2050)iat1mc(i3),iat2mc(i3),vmcha(i3)
      end do
      end if

      ims=0
      do i1=1,na
      if (ims.eq.0.and.imove(i1).eq.0) then
      if1=i1
      ims=1
      end if
      if (ims.eq.1.and.imove(i1).eq.1) then
      write (nunit1,2060)if1,i1-1
      ims=0
      end if
      end do
      if (ims.eq.1) then
      write (nunit1,2060)if1,na
      end if

*     if (qr.eq.'F'.or.qr.eq.'5'.or.qr.eq.'Y') 
      if (ibity.eq.2)
     $write (nunit1,2100)axiss(1),axiss(2),axiss(3),angles(1),
     $angles(2),angles(3)

      if (nrestra.gt.0) write (nunit1,2300)
      do i2=1,nrestra
      write (nunit1,2400)
     $irstra(i2,1),irstra(i2,2),rrstra(i2),
     $vkrstr(i2),vkrst2(i2),rrcha(i2),itstart(i2),itend(i2)
      end do

      if (nrestrav.gt.0) write (nunit1,2500)
      do i2=1,nrestrav
      write (nunit1,2600)
     $irstrav(i2,1),irstrav(i2,2),irstrav(i2,3),
     $vrstra(i2),vkrv(i2),vkr2v(i2),zero
      end do

      if (nrestrat.gt.0) write (nunit1,2700)
      do i2=1,nrestrat
      write (nunit1,2800)
     $irstrat(i2,1),irstrat(i2,2),irstrat(i2,3),
     $irstrat(i2,4),trstra(i2),vkrt(i2),
     $vkr2t(i2),zero
      end do

      if (nrestram.gt.0) write (nunit1,2810)
      do i2=1,nrestram
      write (nunit1,2820)
     $qdir(irstram(i2,1)),irstram(i2,2),irstram(i2,3),
     $rmstra1(i2),irstram(i2,4),irstram(i2,5),rmstra2(i2),
     $rmstra3(i2),rmcha(i2)
      end do

      if (icgeo.eq.0.and.ingeo.eq.0) write (nunit1,2830)
      if (icgeo.eq.1.and.ingeo.eq.1) write (nunit1,2840)
      if (ifreq.eq.1) write (nunit1,2850)
      write (nunit1,2900)
      do i2=1,na
      write (nunit1,3000)i2,qa(i2),c(i2,1),c(i2,2),c(i2,3),
     $qa(i2),irom,irom,chgbgf(i2)
      end do
      write (nunit1,3100)
      if (nsurp.lt.2) then
      do i1=1,na
      write (nunit1,3200)i1,(iag(i1,2+i2),i2=1,iag(i1,2))
      end do
      write (nunit1,3300)
      write (nunit1,3400)estrc
      end if

      write (nunit1,3500)
      write (nunit1,*)
     
      return
 1500 format ('BIOGRF',i4)
 1600 format ('XTLGRF',i4)
 1700 format ('DESCRP ',a60)
 1800 format ('REMARK ',a60)
 1900 format ('FFIELD ',a40)
 2000 format ('RUTYPE NORMAL RUN')
 2010 format ('RUTYPE ENDPO',f6.3)
 2020 format ('RUTYPE MAXIT',i6)
 2030 format ('RUTYPE MAXMOV',i6)
 2033 format ('RUTYPE NO CELL OPT')
 2034 format ('RUTYPE CELL OPT',i6)
 2035 format ('RUTYPE REDO',i6)
 2036 format ('RUTYPE CHARGEMET',i6)
 2040 format ('VCHANGE',f8.4)
 2050 format ('MOLCHARGE',2i4,f6.2)
 2060 format ('FIXATOMS',2i6)    
 2100 format ('CRYSTX ',6f11.5)
 2200 format ('CELLS ',6i5)
 2300 format ('#              At1 At2   R12    Force1  Force2  ',
     $'dR12/dIter(MD) Start (MD) End (MD)')
 2400 format ('BOND RESTRAINT ',2i4,f8.4,f8.2,f8.4,1x,f10.7,2i8)
 2500 format ('#               At1 At2 At3 Angle   Force1  Force2',
     $'  dAngle/dIteration (MD only)')
 2600 format ('ANGLE RESTRAINT ',3i4,2f8.2,f8.4,f9.6)
 2700 format ('#                 At1 At2 At3 At3 Angle   Force1  ',
     $'Force2  dAngle/dIteration (MD only)')
 2800 format ('TORSION RESTRAINT ',4i4,2f8.2,f8.4,f9.6)
 2810 format ('#              x/y/z At1 At2    R   At3 At4 Force1',
     $'  Force2  dR/dIteration (MD only)')
 2820 format ('MASCEN RESTRAINT ',a1,1x,2i4,f8.2,2i4,2f8.2,f9.6)
 2830 format ('GEOUPD')
 2840 format ('NO GEOUPD')
 2850 format ('FREQUENCY')
 2900 format ('FORMAT ATOM   (a6,1x,i5,1x,a5,1x,a3,1x,a1,1x,a5,',
     $'3f10.5,1x,a5,i3,i2,1x,f8.5)')
 3000 format ('HETATM',1x,i5,1x,a2,3x,1x,3x,1x,1x,1x,5x,3f10.5,1x,
     $a5,i3,i2,1x,f8.5)
 3100 format ('FORMAT CONECT (a6,12i6)')
 3200 format ('CONECT',12i6)
 3300 format ('UNIT ENERGY   kcal')
 3400 format ('ENERGY',5x,f14.6)
 3500 format ('END')

      end
**********************************************************************
**********************************************************************

      subroutine writeen(tottime,sum1,sdev,sdeva,sum12,sumt,sump,
     $sumtt,tmax,eaver,eav2,eav3,etot2,ediff)
**********************************************************************   
#include "cbka.blk"
#include "cbkcha.blk"
#include "cbkenergies.blk"
#include "cbkrestr.blk"
#include "cbktorang.blk"
#include "cbktorsion.blk"
#include "cbktregime.blk"
#include "control.blk"
#include "small.blk"

      dimension disres(mrestra)
**********************************************************************
*                                                                    *
*     Write out MD statistics to units 71,73 and 76                  *
*                                                                    *
**********************************************************************
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In writeen'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

      if (nrep1.gt.1)
     $sdev=sqrt((sum12-sum1*sum1/float(nrep1))/float(nrep1-1))
      eavn=eaver/float(mdstep)
      if (mdstep.gt.1)
     $sdeva=sqrt((eav3-eav2*eav2/float(mdstep))/float(mdstep-1))
C      open (71,file='fort.71',status='unknown',access='append')
C      open (73,file='fort.73',status='unknown',access='append')
      write (71,'(i8,2i4,1x,19(f10.2,1x))')mdstep+nprevrun,nmolo,
     $nmolo5,estrc,ekin,estrc+ekin,tempmd,sum1/float(nrep1),eavn,
     $sumt/float(nrep1),tmax,sump/float(nrep1),sdev,sdeva,tset,
     $tstep*1d+15,rmsg,tottime
      write (73,'(i8,1x,14(f10.2,1x))')mdstep+nprevrun,eb,ea,elp,
     $emol,ev,ecoa,ehb,et,eco,ew,ep,ech,efi
      close (71)
      close (73)
 
      if ((sumt/float(nrep1)).gt.tset) then
      if (invt.eq.0) write (*,*)'Switched to NVT in iteration',mdstep
      invt=1
      end if
 
C      if (nrestra.gt.0.or.nrestrat.gt.0)
C     $open (76,file='fort.76',status='unknown',access='append')
 
      if (nrestra.gt.0) then
      do i2=1,nrestra
      call dista2(irstra(i2,1),irstra(i2,2),disres(i2),dx,dy,dz)
      end do
C      open (76,file='fort.76',status='unknown',access='append')
      write (76,'(i8,1x,40f12.4)')mdstep,eres,estrc,
     $(rrstra(i2),disres(i2),i2=1,nrestra)
      end if
 
      if (nrestrat.gt.0) then
C      open (76,file='fort.76',status='unknown',access='append')
      do i2=1,nrestrat
      do i3=1,ntor
      ih1=irstrat(i2,1)
      ih2=irstrat(i2,2)
      ih3=irstrat(i2,3)
      ih4=irstrat(i2,4)
      if (ih1.eq.it(i3,2).and.ih2.eq.it(i3,3).and.ih3.eq.it(i3,4)
     $.and.ih4.eq.it(i3,5)) ittr=i3
      end do
      write (76,'(i8,1x,40f12.4)')mdstep,eres,
     $trstra(i2),thg(ittr)
      end do
      end if
 
      if (nrestra.gt.0.or.nrestrat.gt.0) close(76)
 
      if (nrestram.gt.0) then
C      open (76,file='fort.76',status='unknown',access='append')
      do i2=1,nrestram
      write (76,'(2i8,1x,20f12.4)')mdstep,i2,eres,rmstra1(i2),
     $dismacen(i2)
      end do
      close (76)
      end if

      return
      end
**********************************************************************
************************************************************************

      subroutine molanal

************************************************************************
#include "cbka.blk"
#include "cbkbo.blk"
#include "cbkconst.blk"
#include "cbkdcell.blk"
#include "cbkff.blk"
#include "cbkia.blk"
#include "cbkrbo.blk"
#include "control.blk"
#include "opt.blk"
#include "small.blk"
#include "cbksrtbon1.blk"
      dimension iam(nat,mbond+3),nmolata(nmolmax,nat)
      dimension molfra(nmolmax,nsort),ndup(nmolmax)
      character*40 qmolan1
      character*100 qmolan
      logical found
************************************************************************
*                                                                      *
*     Analyse and output molecular fragments                           *
*                                                                      *
************************************************************************
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In molanal'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

      do i1=1,nmolmax
      do i2=1,nsort
      molfra(i1,i2)=0
      end do
      ndup(i1)=1
      end do

      do i1=1,na
      do i2=1,mbond+3
      iam(i1,i2)=0
      end do
      end do
************************************************************************
*                                                                      *
*     Create connection table based on corrected bond orders           *
*                                                                      *
************************************************************************
      do i1=1,nbon
      if (bo(i1).gt.cutof3) then
      j1=ib(i1,2)
      j2=ib(i1,3)
      iam(j1,2)=iam(j1,2)+1
      iam(j1,2+iam(j1,2))=j2
      iam(j2,2)=iam(j2,2)+1
      iam(j2,2+iam(j2,2))=j1
      end if
      end do
**********************************************************************
*                                                                    *
*     Find molecules                                                 *
*                                                                    *
**********************************************************************
      nmolo6=0
      found=.FALSE.
      DO 61 k1=1,na
      IF (iam(K1,3+mbond).EQ.0) found=.TRUE.
   61 IF (iam(K1,3+mbond).GT.nmolo6) nmolo6=iam(K1,3+mbond)
      IF (.NOT.FOUND) GOTO 62
************************************************************************
*                                                                      *
*     Molecule numbers are assigned. No restrictions are made for the  *
*     sequence of the numbers in the connection table.                 *
*                                                                      *
************************************************************************
      N3=1
   64 N2=N3
      nmolo6=nmolo6+1
      if (nmolo6.gt.nmolmax) stop 'Too many molecules in system'
      iam(N2,3+mbond)=nmolo6
   67 FOUND=.FALSE.
      DO 66 N1=N2+1,na
      IF (iam(N1,3+mbond).NE.0) GOTO 66
      DO 65 L=1,mbond
      IF (iam(N1,l+2).EQ.0) GOTO 66
      IF (iam(iam(N1,l+2),3+mbond).EQ.nmolo6) THEN
      FOUND=.TRUE.
      iam(N1,3+mbond)=nmolo6
      GOTO 66
      ENDIF
   65 CONTINUE
   66 CONTINUE
      IF (FOUND) GOTO 67
      DO 63 N3=N2+1,NA
   63 if (iam(N3,3+mbond).eq.0) goto 64
************************************************************************
*                                                                      *
*     The assigned or input molecule numbers are checked for their     *
*     consistency.                                                     *
*                                                                      *
************************************************************************
   62 FOUND=.FALSE.
      DO 72 N1=1,NA
      DO 71 L=1,mbond
      IF (iam(N1,L+2).EQ.0) GOTO 72
      IF (iam(iam(N1,L+2),3+mbond).NE.iam(N1,3+mbond)) THEN
      FOUND=.TRUE.
      ENDIF
   71 CONTINUE
   72 CONTINUE
      IF (FOUND) THEN
      write (7,'(i4,a40)')na,qmol
      do i1=1,na
      write (7,'(40i4)')i1,iam(i1,1),(iam(i1,2+i2),i2=1,nsbmax),
     $iam(i1,3+mbond)
      end do
      STOP' Mol.nrs. not consistent; maybe wrong cell parameters'
      ENDIF

      do i1=1,nmolo6
      natmol=0
      do i2=1,na
      if (iam(i2,3+mbond).eq.i1) then
      natmol=natmol+1
      nmolata(i1,natmol+1)=i2
      end if
      end do
      nmolata(i1,1)=natmol
      end do
************************************************************************
*                                                                      *
*     Analyze molecules                                                *
*                                                                      *
************************************************************************
      do i1=1,nmolo6
      do i2=1,nmolata(i1,1)
      i3=nmolata(i1,1+i2)
      ityp=ia(i3,1)
      molfra(i1,ityp)=molfra(i1,ityp)+1
      end do
      end do

      do i1=1,nmolo6
      isee=0
      do i2=1,nmolo6
      isee2=1
      do i3=1,nso
      if (molfra(i1,i3).ne.molfra(i2,i3)) isee2=0
      end do
      if (isee2.eq.1.and.i1.gt.i2.and.isee.eq.0) then  !molecule type already exists
      ndup(i2)=ndup(i2)+1
      ndup(i1)=0
      isee=1
      end if

      end do
      end do

C      open (45,file='molfra.out',status='unknown',access='append')
      if (mdstep.eq.0) write (45,100)cutof3
      write (45,110)
      ntotmol=0
      ntotat=0
      vtotmass=zero
      do i1=1,nmolo6
      if (ndup(i1).gt.0) then
*     write (45,110)i1,(molfra(i1,i2),i2=1,nso),ndup(i1)
      ntotmol=ntotmol+ndup(i1)
      qmolan=' '
      qmolan1=' '
      istart=-4
      ihulp=0
      vmass=zero
      do i2=1,nso
      vmass=vmass+molfra(i1,i2)*amas(i2)
      ntotat=ntotat+molfra(i1,i2)*ndup(i1)
      if (molfra(i1,i2).gt.0) then
      istart=istart+6
      iend=istart+5
      if (molfra(i1,i2).gt.1) then
      write (qmolan(istart:iend),'(a2,i3)')qas(i2),molfra(i1,i2)
      else
      write (qmolan(istart:iend-2),'(a2)')qas(i2)
      end if
      end if
      end do
      ihulp=1
      do i2=1,iend
      if (qmolan(i2:i2).ne.' ') then
      qmolan1(ihulp:ihulp)=qmolan(i2:i2)
      ihulp=ihulp+1
      end if
      end do

*     write (45,120)ndup(i1),qmolan(1:iend),vmass
      write (45,120)mdstep,ndup(i1),qmolan1,vmass
      vtotmass=vtotmass+ndup(i1)*vmass
      end if
      end do
      write (45,*)'Total number of molecules:',ntotmol
      write (45,*)'Total number of atoms:',ntotat
      write (45,*)'Total system mass:',vtotmass
      close (45)
      return
  100 format('Bond order cutoff:',f6.4)
  110 format('Iteration Freq. Molecular formula',15x,'Molecular mass')
  120 format(i8,i4,' x  ',a35,f10.4)
      end
************************************************************************
************************************************************************

      subroutine stranal(istart,iend,vout,iout,icheck)

************************************************************************
#include "cbka.blk"
#include "cbkconst.blk"
#include "opt.blk"

      character*1 qchar
      dimension qchar(5)
**********************************************************************
*                                                                    *
*     Analyze string for special characters; find words in string    *
*                                                                    *
**********************************************************************
      qchar(1)=' '
      qchar(2)='/'
      
      ifound1=0
      do i1=istart,200
      ifound2=0
      do i2=1,icheck

      if (qstrana1(i1:i1).eq.qchar(i2)) then
      ifound2=1
      if (ifound1.eq.1) then     !End of word
      iend=i1
      goto 10
      end if

      end if

      end do

      if (ifound2.eq.0.and.ifound1.eq.0) then     !Start of word
      istart2=i1
      ifound1=1
      end if

      end do

   10 continue
      qstrana2=' '
      vout=zero
      iout=0

      if (ifound1.eq.1) then
      qstrana2=qstrana1(istart2:iend-1)
      istart=istart2
      vout=zero
      read (qstrana2,*,end=20,err=20) vout
   20 iout=int(vout)
      end if

      return
      end
************************************************************************
********************************************************************** 

      subroutine dipmom(naold,dpmm,xdip,ydip,zdip,xdir,ydir,zdir)

********************************************************************** 
#include "cbka.blk"
#include "cbkc.blk"
#include "cbkch.blk"
#include "cbkconst.blk"
#include "control.blk"
#include "small.blk"
********************************************************************** 
*                                                                    *
*     Calculate and output dipole moment                             *
*                                                                    *
********************************************************************** 
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In dipmom'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if
************************************************************************
*                                                                      *
*     CONVERSION FACTOR TO DEBYE UNITS IS CALCULATED                   *
*     THE CALCULATION IS INITIALIZED                                   *
*                                                                      *
************************************************************************
 
      ELCHG=1.60217733D-19       ! [C]       = [As]
      CLIGHT=2.99792458D8        !             [m/s]
      DBCONV=ONE/(CLIGHT*ELCHG*1.0D11)
 
      CHCPX=ZERO
      CHCPY=ZERO
      CHCPZ=ZERO
      CHCMX=ZERO
      CHCMY=ZERO
      CHCMZ=ZERO
      XDIP=ZERO
      YDIP=ZERO
      ZDIP=ZERO
      XGRD=ZERO
      YGRD=ZERO
      ZGRD=ZERO
************************************************************************
*                                                                      *
*     CALCULATION OF MAGNITUDE AND CENTRES OF POSITIVE AND NEGATIVE    *
*     CHARGES                                                          *
*                                                                      *
************************************************************************
 
      if (na.eq.0) na=naold
      CHRG=ZERO
      DO 4 K1=1,NA
      CHK1=CH(K1)
      IF (CHK1.EQ.ZERO) GOTO 4
      IF (CHK1.LT.ZERO) GOTO 3
      CHRG=CHRG+CHK1
      CHCPX=CHCPX+CHK1*C(K1,1)
      CHCPY=CHCPY+CHK1*C(K1,2)
      CHCPZ=CHCPZ+CHK1*C(K1,3)
      GOTO 4
    3 CHCMX=CHCMX-CHK1*C(K1,1)
      CHCMY=CHCMY-CHK1*C(K1,2)
      CHCMZ=CHCMZ-CHK1*C(K1,3)
    4 CONTINUE
 
************************************************************************
*                                                                      *
*     CALCULATION OF DISTANCE BETWEEN CENTRES AND OF DIPOLE MOMENT     *
*     IN DEBIJE UNITS                                                  *
*                                                                      *
************************************************************************
 
      CHDSTX=CHCPX-CHCMX
      CHDSTY=CHCPY-CHCMY
      CHDSTZ=CHCPZ-CHCMZ
      DPMM=SQRT(CHDSTX*CHDSTX+CHDSTY*CHDSTY+CHDSTZ*CHDSTZ)/DBCONV
      IF(DPMM.LT.1.0D-4)RETURN
      XDIP=HALF*(CHCPX+CHCMX)/CHRG
      YDIP=HALF*(CHCPY+CHCMY)/CHRG
      ZDIP=HALF*(CHCPZ+CHCMZ)/CHRG
      GRTST=MAX(CHDSTX,CHDSTY,CHDSTZ)
      XDIR=-CHDSTX/GRTST
      YDIR=-CHDSTY/GRTST
      ZDIR=-CHDSTZ/GRTST
      open (64,file='dipole.out',status='unknown')
      write (64,100)dpmm,xdip,ydip,zdip,xdir,ydir,zdir
      close (64)
 
  100 format ('Dipole moment (Debye):',f12.4,' Location:',3f12.4,
     $' Direction (-side):',3f12.4) 
      return
      end
************************************************************************
**********************************************************************

      subroutine readtraj(ivels)

**********************************************************************
#include "cbka.blk"
#include "cbkatomcoord.blk"
#include "cbkc.blk"
#include "cbkconst.blk"
#include "cbkdistan.blk"
#include "cbktregime.blk"
#include "cellcoord.blk"
#include "control.blk"
#include "small.blk"
#include "cbkinit.blk"
**********************************************************************
*                                                                    *
*     Read in trajectory file                                        *
*                                                                    *
**********************************************************************
c$$$      if (ndebug.eq.1) then
c$$$C      open (65,file='fort.65',status='unknown',access='append')
c$$$      write (65,*) 'In readtraj'
c$$$      call timer(65)
c$$$      close (65)
c$$$      end if

      open(unit=66,file='moldyn.vel',status='old',err=10)
      ivels=1
      read (66,*)
      read (66,100)aaxis,baxis,caxis
      read (66,100)angles(1),angles(2),angles(3)
      if (qr.eq.'F'.or.qr.eq.'P'.or.ngeofor.eq.1) then
      axis(1)=aaxis
      axis(2)=baxis
      axis(3)=caxis
      axiss(1)=axis(1)
      axiss(2)=axis(2)
      axiss(3)=axis(3)
      angle(1)=angles(1)
      angle(2)=angles(2)
      angle(3)=angles(3)
      halfa=angle(1)*dgrrdn
      hbeta=angle(2)*dgrrdn
      hgamma=angle(3)*dgrrdn
      sinalf=sin(halfa)
      cosalf=cos(halfa)
      sinbet=sin(hbeta)
      cosbet=cos(hbeta)
      cosphi=(cos(hgamma)-cosalf*cosbet)/(sinalf*sinbet)
      if (cosphi.gt.1.0) cosphi=1.0
      sinphi=sqrt(one-cosphi*cosphi)
      tm11=axis(1)*sinbet*sinphi
      tm21=axis(1)*sinbet*cosphi
      tm31=axis(1)*cosbet
      tm22=axis(2)*sinalf
      tm32=axis(2)*cosalf
      tm33=axis(3)
      end if
      if (aaxis.ne.axis(1).or.baxis.ne.axis(2).or.caxis.ne.axis(3))
     $stop 'Wrong cell parameters in moldyn.vel'
      read (66,200)nan
      if (nan.ne.na) stop 'Wrong number of atoms in moldyn.vel-file'
      if (nbiolab.eq.1) write (*,*)'Warning: using labels in vels-file'
      read (66,250)((c(i,j),j=1,3),qlabel(i),i=1,na)
      read (66,*)
      read (66,300)((vel(j,i),j=1,3),i=1,na)
      read (66,*)
      read (66,300)((accel(j,i),j=1,3),i=1,na)
      read (66,*)
      read (66,300,end=10,err=10)((aold(j,i),j=1,3),i=1,na)
      read (66,*)
      read (66,300,end=10,err=10)tempmd
      read (66,*)
      read (66,350,end=10,err=10)nsbma2
   10 continue
**********************************************************************
*                                                                    *
*     Format part                                                    *
*                                                                    *
**********************************************************************
  100 format(3d15.8)
  200 format(i4)
  250 format(3d24.15,1x,a5)
  300 format(3d24.15)
  350 format(i3)
  400 format (8i3,8f8.4)
      return
      end
**********************************************************************