CHARMM Element source/gamint/gamini.src 1.1
      SUBROUTINE GAMINI(COMLYN,COMLEN)
C-----------------------------------------------------------------------
C     Define the set of quantum mechanical atoms and set up the data
C     structures.
C
C     April 1993,   Milan Hodoscek,  Pure QM part
C     October 1993,   -- " --        QM/MM interface
C
C     Currently there is no cutoff method incorporated, but this should
C     be straightforward. No big loss since Natom(QM) is usually small so 
C     Natom(QM)*Natom(MM)/2 is also small. Also there shouldn't be any 
C     inconsistency with MM cutoff methods, since these are different
C     terms anyway. 
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
##INCLUDE '~/charmm_fcm/code.fcm'
##INCLUDE '~/charmm_fcm/coord.fcm'
##INCLUDE '~/charmm_fcm/psf.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/energy.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      character*240   name
      integer kstrm,koutu
      COMMON /CHGMIO/  KSTRM,KOUTU,name
C
      CHARACTER COMLYN*(*)
      INTEGER   COMLEN
C
      INTEGER   ISLCT,OLDUSD
      LOGICAL   QQINP
CCCC  test code...
      logical mopac
      common /mpctst/mopac
C
##IFN CADPAC
##IFN GAMESS
      CALL WRNDIE(-1,'<GAMINI>','GAMESS code not compiled.')
##ELSE
C
      OLDUSD = LSTUSD
      ISLCT  = ALLSTK(INTEG4(NATOM))
      CALL SELCTA(COMLYN,COMLEN,STACK(ISLCT),X,Y,Z,WMAIN,.TRUE.)
      QGMREM=(INDXA(COMLYN,COMLEN,'REMO').GT.0)
      IF(PRNLEV.GE.2) THEN
         IF(QGMREM) THEN
            WRITE(OUTU,22)
     &        'REMOve: Classical energies within QM atoms are removed.'
         ELSE
            WRITE(OUTU,22)
     &    'No REMOve: Classical energies within QM atoms are retained.'
         ENDIF
      ENDIF
C
      QGMEXG=(INDXA(COMLYN,COMLEN,'EXGR').GT.0)
      IF(PRNLEV.GE.2) THEN
         IF(QGMEXG) THEN
            WRITE(OUTU,22)
     &     'EXGRoup: QM/MM Electrostatics for link host groups removed.'
         ELSE
            WRITE(OUTU,22)
     &     'No EXGRoup: QM/MM Elec. for link atom host only is removed.'
         ENDIF
      ENDIF
C
      QBLUCH=(INDXA(COMLYN,COMLEN,'BLUR').GT.0)
      IF(QBLUCH) THEN
         IF(PRNLEV.GE.2) THEN
            WRITE(OUTU,22)
     &      'BLUR: Blurred charges will be used on some atoms.'
         ENDIF
      ENDIF
C
      QNOGU=(INDXA(COMLYN,COMLEN,'NOGU').GT.0)
      IF(QNOGU) THEN
         KGUES=0
         IF(PRNLEV.GE.2) THEN
            WRITE(OUTU,22)
     &      'NOGUess: Initial guess obtained from previous step.'
         ENDIF
      ENDIF
C
      QQINP=(INDXA(COMLYN,COMLEN,'QINP').GT.0)
      IF(PRNLEV.GE.2) THEN
         IF(QQINP) THEN
            WRITE(OUTU,22)
     &     'QINP: Charges will input for QM atoms.'
         ELSE
            WRITE(OUTU,22)
     &     'No QINP: Charges will be based on atomic numbers.'
         ENDIF
      ENDIF
      MOPAC=(INDXA(COMLYN,COMLEN,'MOPAC').GT.0)
      IF(PRNLEV.GE.2) THEN
         IF(MOPAC) THEN
            WRITE(OUTU,22)
     &     'MOPAC: A MOPAC calculation will be performed.'
         ENDIF
      ENDIF
  22  FORMAT('GAMINT> ',A)
C
      MUSTUP=.TRUE.
C
      CALL COPSEL(STACK(ISLCT),QQINP)
C
      CALL CH2GMS(.TRUE.)
C
C     This initialize gamess data 
      QINIGM=.TRUE.
      CALL GAMESS
C
C     Report QM/MM repulsion energy, also when no derivatives involved
      IF(PRNLEV.GE.2) CALL CGREPE(NATOM)
C
C     Now we have to put pointer for input file at the end of gamess input
CCC      CALL GAMEND(99)
C
      CALL FRESTK(LSTUSD - OLDUSD)
      COMLEN = 0
C
      RETURN
      END
C
      SUBROUTINE CH2GMS(BFIRST)
C-----------------------------------------------------------------------
C     Define CHARMM atoms as point charges and copy to COMMON/CHMGMS/
C
C     To simplify changes in GAMESS we build contiguous arrays:
C     (XCHM,YCHM,ZCHM,QCHM) for .not.QM atoms.
C     [NOTE: Similar code is present also in CGREP for
C            nuclear repulsion derivatives]
C
C     On QM atoms we don't care for cutoff, we just take all of the
C     atoms in the system. This is not inconsistent with MM since QMs
C     are special anyway. Also calculation time is linear with number
C     of MM atoms!
C
C     [NOTE: It should be straightforward to implement usage of
C            variety of cutoff methods implemented in the CHARMM:
C            Just use nonbond array here and disable CGREP routine;
C            but you need to specify charges on QM atoms as their
C            atomic number in the RTF file!!!
C
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
##INCLUDE '~/charmm_fcm/coord.fcm'
##INCLUDE '~/charmm_fcm/psf.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
##INCLUDE '~/charmm_fcm/consta.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
C
      LOGICAL BFIRST
C
      INTEGER I,N,NBLUR
      REAL*8 TPIPOH,SIGM1,SIGM2,SIGM3
      REAL*8 RBR
C
C
C     All atoms which are not quantum contribute to electrostatic
C     interaction in the QM part. See MJF reference:
C     J. Comp. Chem., Vol. 11, No. 6, 700-733 (1990)
C
      N=0
      NBLUR=0
C     tpipoh = 2*sqrt(PI**3)
##IF SINGLE
      TPIPOH=11.136655993663416E0
##ELSE
      TPIPOH=11.136655993663416D0
##ENDIF
      RBR=ONE/BOHRR
      CALL FILLI4(IBLUCH,NATOM,0)
      DO I = 1,NATOM
         IF(QBLUCH) THEN
            IF (IGMSEL(I) .LE. 0) THEN
               N = N + 1
               XCHM(N)=X(I)*RBR
               YCHM(N)=Y(I)*RBR
               ZCHM(N)=Z(I)*RBR
               QCHM(N)=CG(I)
C
C     Fill in blurred charges arrays
C
               IF (ABS(WMAIN(I)).GE.RSMALL) THEN
                  IF(WMAIN(I).GT.NINE99) THEN
                     QCHM(N)=ZERO
                  ELSE
                     NBLUR=NBLUR+1
                     IBLUCH(NBLUR)=N
                     SIGM1=BOHRR/WMAIN(I)
                     SIGM2=SIGM1*SIGM1
                     SIGM3=SIGM2*SIGM1
                     EBLUCH(NBLUR)=SIGM2
                     IF(BFIRST) THEN
                        CGBLCH(NBLUR)=CG(I)
                     ENDIF
                     SGBLCH(NBLUR)=WMAIN(I)*RBR
CCC                     CBLUCH(NBLUR)=CGBLCH(NBLUR)*SIGM3/TPIPOH
                     CBLUCH(NBLUR)=CGBLCH(NBLUR)*SIGM3*0.5079490874739
C
C                    Put to zero both MM charges
C                    a) the one which goes to GAMESS (QCHM)
C                    b) the one which goes to CHARMM (CG)
C
                     CG(I)=ZERO
                     QCHM(N)=ZERO
C
C                    For now put back MM charges
C                    They are dealt in CHARMM, and we only have
C                    Blurred charges interacting with QM atoms.
C                    No blur-blur! This is done in standard CHARMM way!
C
                     CG(I)=CGBLCH(NBLUR)
                  ENDIF
               ENDIF
            ENDIF
         ELSE
            IF (IGMSEL(I) .EQ. 0) THEN
               N = N + 1
               XCHM(N)=X(I)/BOHRR
               YCHM(N)=Y(I)/BOHRR
               ZCHM(N)=Z(I)/BOHRR
               QCHM(N)=CG(I)
CCCC           write(*,*)'ch2gms>',N,I,XCHM(N),YCHM(N),ZCHM(N),QCHM(N)

            ENDIF
         ENDIF
      ENDDO
C
      NCHMAT=N
      NBLUCH=NBLUR
C
      RETURN
      END
C
      SUBROUTINE CHMDAT(AATOM,AZNUC,CORD,NAT)
C-----------------------------------------------------------------------
C     Define the set of quantum mechanical atoms and set up the data
C     structures.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
##INCLUDE '~/charmm_fcm/coord.fcm'
##INCLUDE '~/charmm_fcm/psf.fcm'
##INCLUDE '~/charmm_fcm/rtf.fcm'
##INCLUDE '~/charmm_fcm/param.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
C
C     COMMON from GAMESS
C
C     The following is defined in GAMESS after call to CHMDAT
C      REAL*8 HMLTN
C      LOGICAL MPCWFN
C      COMMON /MPCLNK/ HMLTN,MPCWFN
CCCC  test code...
      logical mopac
      common /mpctst/mopac
C
      CHARACTER*10 AATOM(MAXGMS)
      REAL*8 AZNUC(MAXGMS), CORD(MAXGMS,3)
      INTEGER NAT
C
      INTEGER I,N,NSLCT,NATMM,NATQM,NATLNK
      CHARACTER*6 ELE
C
      N=0
      DO I = 1,NATOM
        IF ((IGMSEL(I) .GE. 1) .OR. MOPAC ) THEN
           N = N + 1
           CORD(N,1)=X(I)
           CORD(N,2)=Y(I)
           CORD(N,3)=Z(I)
           CUNIQ(N,1)=CORD(N,1)
           CUNIQ(N,2)=CORD(N,2)
           CUNIQ(N,3)=CORD(N,3)
C
C          Added this condition for lone pairs.
C          This prevents to change QM region within one run
C
           IF (QINIGM) THEN
              AATOM(N)=TYPE(I)
              IF(TYPE(I)(1:3).EQ.'QQH')AATOM(N)='H         '
C
              CALL FINDEL(ATCT(IAC(I)),AMASS(I),I,ELE,AZNUC(N),QINIGM)
C
              IF (IGMSEL(I) .LT. 1 ) THEN
C This is a sparkle
                 AZNUC(N)=1
                 AATOM(N)='+         '
                 write(*,*)'CHMDAT>igmsel,n,i=',igmsel(i),n,i
              ENDIF
              UATOM(N)=AATOM(N)
              UZNUC(N)=AZNUC(N)
CCC           write(*,*)'CHMDAT>',i,n,type(i),aatom(n),aznuc(n)
           ENDIF
        ENDIF
      ENDDO
C
      NAT=N
      NATREL=NAT
      UATOM(NAT+1)='$END      '
C
      IF (NAT .LE. 0) CALL WRNDIE(0,'<CHMDAT>',
     $                   'No quantum mechanical atoms selected.')
      NATMM = NATOM - N
      NATQM = N
      NGAMES = N
C
      NSLCT = 0
      DO I = 1,NATOM
        IF (IGMSEL(I).EQ.2) NSLCT = NSLCT + 1
      ENDDO
      NATLNK = NSLCT
C
C     Write out some information and options requested.
C
      IF(QBLUCH.AND.PRNLEV.GE.2.AND.QINIGM) WRITE (OUTU,'(/,8X,A,I5)')
     $  ' The number of blurred MM charges         = ',NBLUCH
      IF ((PRNLEV.GE.2).AND.QINIGM) WRITE (OUTU,'(/,1X,A,/)')
     $  ' GAMDFN> Some atoms will be treated quantum mechanically.'
      IF ((PRNLEV.GE.2).AND.QINIGM) WRITE (OUTU,'(4(8X,A,I5,/),/)')
     $  ' The number of quantum   mechanical atoms = ',NATQM,
     $  ' The number of molecular mechanical atoms = ',NATMM,
     $  ' The number of MM atoms excluded from QM  = ',NATMM-NCHMAT,
     $  ' Of which the number of QM/MM link atoms  = ',NATLNK
C
CCC      write(*,*)'x',(Cord(i,1),i=1,natom)
CCC      write(*,*)'y',(Cord(i,2),i=1,natom)
CCC      write(*,*)'z',(Cord(i,3),i=1,natom)
C
      RETURN
      END
C
      SUBROUTINE GMSENE(GTOT,X,Y,Z,DX,DY,DZ,NATOM)
C-----------------------------------------------------------------------
C
C     Get energy and forces from GAMESS
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
      REAL*8 GTOT,X(*),Y(*),Z(*),DX(*),DY(*),DZ(*)
      INTEGER NATOM
C
##INCLUDE '~/charmm_fcm/consta.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
C
      INTEGER ICHARM
C
      REAL*8 E, EG
      COMMON /FUNCT/ E, EG(3*MAXGMS)
C
      INTEGER NAT,ICH,MUL,NUM,NX,NE,NA,NB
      REAL*8 ZAN, C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NX,NE,NA,NB,
     $     ZAN(MAXGMS),C(3,MAXGMS)
C
      REAL*8 HMLTN
      LOGICAL MPCWFN
      COMMON /MPCLNK/ HMLTN,MPCWFN
C
      REAL*8 QMMMRP, RBR
      INTEGER IPT,I,N
      logical mopac
C
C     Are there any QM atoms?
C
      IF(NGAMES.EQ.0) RETURN
C
      RBR=ONE/BOHRR
C Update coordinates
      N=0
      DO I = 1,NATOM
        IF (IGMSEL(I) .GE. 1) THEN
           N = N + 1
           C(1,N)=X(I)*RBR
           C(2,N)=Y(I)*RBR
           C(3,N)=Z(I)*RBR
        ENDIF
      ENDDO
C
      CALL CH2GMS(.FALSE.)
C
      CALL CGREP(NATOM,QMMMRP,DX,DY,DZ)
      IF(PRNLEV.GE.6) WRITE(OUTU,'(A,2F17.6)')
     $     'QM/MM repulsion is (au,kcal) ',QMMMRP,QMMMRP*TOKCAL
C
C     From gamess/grd1.src in case of other WFns.
C
      IF(NCHMAT.NE.0) THEN
         DO ICHARM=1,NCHMAT
            DXELMM(ICHARM)=ZERO
            DYELMM(ICHARM)=ZERO
            DZELMM(ICHARM)=ZERO
         ENDDO
      ENDIF

CCC      do i=1,3*ngames
CCC         eg(i)=zero
CCC      enddo
C
CC      CALL START
CC      irest=0
CC      CALL GRADX
      QINIGM=.FALSE.
      CALL GAMESS
CCC      call gamend(99)
C
C     MOPAC energies are already in kcal/mole.
C
C     This is NOT TRUE ANYMORE!!!!!
C
C      IF(MPCWFN) THEN
C         GTOT=E+QMMMRP*TOKCAL
C      ELSE
C         GTOT=(E+QMMMRP)*TOKCAL
C      ENDIF
C
##IF PARALLEL
      IF(MYNOD.EQ.0) THEN
         GTOT=(E+QMMMRP)*TOKCAL
      ELSE
         GTOT=ZERO
      ENDIF
##ELSE
      GTOT=(E+QMMMRP)*TOKCAL
##ENDIF
C
C     MM atoms, without igmsel(i)=-1 !!
C
      N = 0
      DO I = 1,NATOM
         IF ((IGMSEL(I).EQ.0).OR.((IGMSEL(I).LE.0).AND.QBLUCH)) THEN
            N = N + 1
            DX(I) = DX(I) - DXELMM(N)*TOKCAL*RBR
            DY(I) = DY(I) - DYELMM(N)*TOKCAL*RBR
            DZ(I) = DZ(I) - DZELMM(N)*TOKCAL*RBR
C
            IF(PRNLEV.GT.6) THEN
               WRITE(OUTU,334) I,N, - DXELMM(N)*TOKCAL*RBR
     &              , - DYELMM(N)*TOKCAL*RBR
     &              , - DZELMM(N)*TOKCAL*RBR
            ENDIF
C
         ENDIF
      ENDDO
C
C
C     QM atoms
C
##IF PARALLEL
      IF(MYNOD.GT.0) RETURN
##ENDIF
      N = 0
      DO I = 1,NATOM
         IF (IGMSEL(I) .GE. 1) THEN
            N = N + 1
            IPT=3*(N-1)+1
            DX(I) = DX(I) + EG(IPT)*TOKCAL*RBR
            DY(I) = DY(I) + EG(IPT+1)*TOKCAL*RBR
            DZ(I) = DZ(I) + EG(IPT+2)*TOKCAL*RBR
C
            IF(PRNLEV.GT.6) THEN
               WRITE(OUTU,334) I,N, EG(IPT)*TOKCAL*RBR    
     &              , EG(IPT+1)*TOKCAL*RBR
     &              , EG(IPT+2)*TOKCAL*RBR
 334           FORMAT(2I5,3F14.6,3F14.9)
            ENDIF
C
         ENDIF
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE CGREP(NATOM,REPULS,DX,DY,DZ)
C-----------------------------------------------------------------------
C
C     This routine calculates nuclear repulsion between
C     CHARMM atoms and GAMESS atoms + derivatives
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/coord.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
##INCLUDE '~/charmm_fcm/consta.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
C
      REAL*8 REPULS,DX(*),DY(*),DZ(*)
      INTEGER NATOM
C
      INTEGER NAT,ICH,MUL,NUM,NX,NE,NA,NB
      REAL*8 ZAN, C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NX,NE,NA,NB,
     $     ZAN(MAXGMS),C(3,MAXGMS)
C
      INTEGER I,J,N,L,KBLUCH,KKBLCH,NN
      REAL*8 ERRF,ETMP,SIGM1,SIGM2,TSQP
      REAL*8 Q1,Q2,X1,X2,Y1,Y2,Z1,Z2,R12,RR12,EL,ELR
      REAL*8 X12,Y12,Z12,RBR
C
C
C     Put derivatives in the right places in DX,DY,DZ arrays
C     and take Newton's 3rd law into account. Also 
C     here we don't care for cutoff for the same reason as
C     above.
C
      RBR=ONE/BOHRR
      L = 0
      TSQP = TWO/SQRT(PI)
      REPULS = ZERO
C DIRTY ...
      CALL ERRT
C DIRTY ...
      IF(QBLUCH) CALL ERRT
##IF PARALLEL
      IF (MYNOD.GT.0) RETURN
##ENDIF
C
C     This loop is for QM nuclei - MM atoms electrostatic interaction
C     It deals also with QM nuclei - Blurred charge interaction
C
      DO I = 1,NATOM
         IF (IGMSEL(I) .GE. 1) THEN
C
C        QM atoms
C
            L=L+1
            Q1 = ZAN(L)
            IF(FQQCHG(L).GT.-NINE99) Q1 = FQQCHG(L)
            X1 = C(1,L)
            Y1 = C(2,L)
            Z1 = C(3,L)
            N = 0
CCC            write(*,'(''CGREP>x1,y1,z1,q1='',4f12.5)')x1,y1,z1,q1
            KBLUCH=1
            DO J = 1, NATOM
              IF((IGMSEL(J).EQ.0).OR.((IGMSEL(J).LE.0).AND.QBLUCH))THEN
C
C           MM atoms
C
                  N = N + 1
                  X2 = XCHM(N)
                  Y2 = YCHM(N)
                  Z2 = ZCHM(N)
                  X12 = X1-X2
                  Y12 = Y1-Y2
                  Z12 = Z1-Z2
                  RR12 = X12**2+Y12**2+Z12**2
                  R12 = SQRT(RR12)
                  IF(QBLUCH.AND.(N.EQ.IBLUCH(KBLUCH))) THEN
C
C                    QM nuclei - Blurred charge interaction
C
CCC                     write(*,'(''CGREP>r,kbluch,ibluch,cbluch,ebluch='',
CCC     $                    f12.5,3i5,2f12.5)')r12,kbluch,ibluch(kbluch),
CCC     $                    cbluch(kbluch),ebluch(kbluch),sgblch(kbluch)
                     Q2=CGBLCH(KBLUCH)
CCCC                     Q2=CBLUCH(KBLUCH)
                     ETMP=R12/SGBLCH(KBLUCH)
CCC                     write(*,'(''CGREP>s,r,r/s,q2='',4f12.5)')
CCC     $                    SGBLCH(KBLUCH),r12,etmp,q2
CCC                     write(*,'(''CGREP>erf='',f12.5)')errf(etmp)
                     EL = Q1*Q2/R12
                     ELR=EL/RR12
     $                  *(ERRF(ETMP)-TSQP*ETMP*EXP(-ETMP*ETMP))
     $                  *TOKCAL/BOHRR
                     EL = EL*ERRF(ETMP)
                     KBLUCH=KBLUCH+1
                  ELSE
                     Q2 = QCHM(N)
                     EL = Q1*Q2/R12
                     ELR = EL/RR12*TOKCAL/BOHRR
                  ENDIF
                  REPULS = REPULS + EL
CCC                  write(*,'(''CGREP>x2,y2,z2,q2='',4f12.5)')x2,y2,z2,q2
CCC                  write(*,'(''CGREP>repuls='',f20.10)')repuls
##IF FLUCQ
                  CALL FQQCOR(J,EL)
##ENDIF
                  DX(I) = DX(I) - X12*ELR
                  DX(J) = DX(J) + X12*ELR
                  DY(I) = DY(I) - Y12*ELR
                  DY(J) = DY(J) + Y12*ELR
                  DZ(I) = DZ(I) - Z12*ELR
                  DZ(J) = DZ(J) + Z12*ELR
              ENDIF
            ENDDO
         ENDIF
      ENDDO
CCC      write(*,'(''CGREP>repuls='',f20.10)')repuls
C
C     This loop is for Blurred charges - MM atoms electrostatic interaction
C     It deals also with the Blurred charge - Blurred charge interaction
C
C============================
C
C     SIMPLIFICATION: (????)
C     For now we deal with blurred charges not interacting
C     with QM region in the CHARMM as classical atoms.
C     The following code is not good for MM - Blur and Blur - Blur,
C     because it doesn't deal correctly with the bonded atoms!!
C
C     USE: ~/charmm_fcm/inbnd.fcm' '~/charmm_fcm/heap.fcm'
C     USE: heap(bnbnd(inblo)),heap(bnbnd(jnb)))
C
C
      IF(QBLUCH) RETURN
      N = 0
      KBLUCH = 1
      DO I = 1,NATOM
         IF (QBLUCH.AND.(IGMSEL(I) .LE. 0)) THEN
            N = N + 1
            IF(N.EQ.IBLUCH(KBLUCH)) THEN
C
C              Found KBLUCH-th blurred charge.
C
               X1 = XCHM(N)
               Y1 = YCHM(N)
               Z1 = ZCHM(N)
               Q1 = CGBLCH(KBLUCH)
               SIGM1 = SGBLCH(KBLUCH)
               KBLUCH = KBLUCH + 1
               NN = 0
               KKBLCH = 1
               DO J = 1,NATOM
                  IF (IGMSEL(J).LE.0) THEN
C
C                 Either MM or Blurred atom
C
                     NN = NN + 1
                     X2 = XCHM(NN)
                     Y2 = YCHM(NN)
                     Z2 = ZCHM(NN)
                     X12 = X1-X2
                     Y12 = Y1-Y2
                     Z12 = Z1-Z2
                     RR12 = X12**2+Y12**2+Z12**2
                     R12 = SQRT(RR12)
                     ELR=ZERO
                     EL=ZERO
                     IF (NN.EQ.IBLUCH(KKBLCH)) THEN
                        IF (KKBLCH.GE.KBLUCH) THEN
C
C                          This is for blurred charge - blurred charge
C                          Blurred charges with index less than KBLUCH
C                          are ignored (EL,ELR is put to zero)
C
                           SIGM2 = ONE/SQRT(EBLUCH(KKBLCH))
                           ETMP=R12/SQRT(SIGM1*SIGM1+SIGM2*SIGM2)
                           EL=Q1*CGBLCH(KKBLCH)/R12
                           ELR=EL/R12
     $                        *(ERRF(ETMP)-TSQP*ETMP*EXP(-ETMP*ETMP))
     $                        *TOKCAL/BOHRR
                           EL=EL*ERRF(ETMP)
                        ENDIF
                        KKBLCH = KKBLCH + 1
                     ELSE
C                       This is for Blurred charge - MM charge
                        ETMP=R12/SIGM1
                        EL=Q1*QCHM(NN)/R12
                        ELR=EL/R12
     $                     *(ERRF(ETMP)-TSQP*ETMP*EXP(-ETMP*ETMP))
     $                     *TOKCAL/BOHRR
                        EL=EL*ERRF(ETMP)
                     ENDIF
                     REPULS=REPULS+EL
##IF FLUCQ
                     CALL FQQCOR(J,EL)
##ENDIF
                     DX(I) = DX(I) - X12*ELR
                     DX(J) = DX(J) + X12*ELR
                     DY(I) = DY(I) - Y12*ELR
                     DY(J) = DY(J) + Y12*ELR
                     DZ(I) = DZ(I) - Z12*ELR
                     DZ(J) = DZ(J) + Z12*ELR
                  ENDIF
               ENDDO
            ENDIF
         ENDIF
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE CGREPE(NATOM)
C-----------------------------------------------------------------------
C
C     This routine calculates nuclear repulsion between
C     CHARMM atoms and GAMESS atoms - IT DOES NOT WORK WITH BLUR!!!
C                                     Try to get rid of it
C                                     by filling ETERM array with the
C                                     terms calculated in the CGREP
C                                     subroutine
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
##INCLUDE '~/charmm_fcm/consta.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
C
C
      INTEGER NATOM
C
      INTEGER NAT,ICH,MUL,NUM,NX,NE,NA,NB
      REAL*8 ZAN, C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NX,NE,NA,NB,
     $     ZAN(MAXGMS),C(3,MAXGMS)
C
      REAL*8 E, EG
      COMMON /FUNCT/ E, EG(3*MAXGMS)
C
      INTEGER I,J,N,L,KBLUCH
CCC      REAL*8 ERRF
      REAL*8 Q1,Q2,X1,X2,Y1,Y2,Z1,Z2,R12,RR12,EL
      REAL*8 X12,Y12,Z12
      REAL*8 REPULS
C
      L = 0
      REPULS = ZERO
      DO I = 1,NATOM
         IF (IGMSEL(I) .GE. 1) THEN
C
C        QM atoms
C
            L=L+1
            Q1 = ZAN(L)
            IF(FQQCHG(L).GT.-NINE99) Q1 = FQQCHG(L)
            X1 = C(1,L)
            Y1 = C(2,L)
            Z1 = C(3,L)
            N = 0
            KBLUCH=1
            DO J = 1, NATOM
              IF((IGMSEL(J).EQ.0).OR.((IGMSEL(J).LE.0).AND.QBLUCH))THEN
C
C        MM atoms
C
                  N = N + 1
                  X2 = XCHM(N)
                  Y2 = YCHM(N)
                  Z2 = ZCHM(N)
                  X12 = X1-X2
                  Y12 = Y1-Y2
                  Z12 = Z1-Z2
                  RR12 = X12**2+Y12**2+Z12**2
                  R12 = SQRT(RR12)
CCC                  IF(QBLUCH.AND.(N.EQ.IBLUCH(KBLUCH))) THEN
CCC                     Q2=CBLUCH(KBLUCH)*ERRF(R12/SQRT(EBLUCH(KBLUCH)))
CCC                     KBLUCH=KBLUCH+1
CCC                  ELSE
                     Q2 = QCHM(N)
CCC                  ENDIF
                  EL = Q1*Q2/R12
                  REPULS = REPULS + EL
               ENDIF
            ENDDO
         ENDIF
      ENDDO
C
      IF(PRNLEV.GE.2) THEN
         WRITE(OUTU,'(A,2F20.8)')'QM/MM repulsion (a.u.,kcal/mole) = ',
     $        REPULS,REPULS*TOKCAL
         WRITE(OUTU,'(A,2F20.8)')'QM/MM total en. (a.u.,kcal/mole) = ',
     $        E+REPULS,(E+REPULS)*TOKCAL
      ENDIF
C
      RETURN
      END
C
      SUBROUTINE CHGMIU(IR,IW,inize)
C-----------------------------------------------------------------------
C
C     Routine to define unit numbers for input/output of GAMESS
C     The rest of files are opened by gamess with names defined
C     by environment variables set by shell or ENVIron CHARMM command.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
C
      character*240   name
      integer kstrm,koutu
      COMMON /CHGMIO/  KSTRM,KOUTU,name
      logical         qopen,qform,qwrite
C
      INTEGER IR,IW,inize,ip
C
C      The following is the complete table of I/O
C      used in GAMESS program
C      Wich ones are used depend on the type of GAMESS job.
C
C      IRCDATA $SCR/$JOB.irc
C        INPUT $SCR/$JOB.F05
C        PUNCH $SCR/$JOB.dat
C      INTGRLS $SCR/$JOB.F08
C       ORDINT $SCR/$JOB.F09
C       JKFILE $SCR/$JOB.F09
C      DICTNRY $SCR/$JOB.F10  ! ENVI must be always present
C      DRTFILE $SCR/$JOB.F11
C      CIVECTR $SCR/$JOB.F12
C      NTNFMLA $SCR/$JOB.F13
C       CIINTS $SCR/$JOB.F14
C       WORK15 $SCR/$JOB.F15  ! ENVI must be always present
C       WORK16 $SCR/$JOB.F16
C      CSFSAVE $SCR/$JOB.F17
C      FOCKDER $SCR/$JOB.F18
C       DASORT $SCR/$JOB.F20  ! ENVI must be always present
C
      IR=ISTRM
      IW=OUTU
      ip=7
C
      KSTRM=IR
      KOUTU=IW
      if (inize.ne.-1) return
      if(ir.eq.5)ir=99
      if(inize.ge.0) then
         write(*,*)'CHGMIU-0>ir,ip',ir,ip
         CALL SEQOPN(IR,'INPUT', 'OLD',.TRUE., 'FORMATTED')
         CALL SEQOPN(IP,'PUNCH', 'NEW',.FALSE.,'FORMATTED')
CCC   call vinqre('unit',name,240,length,qopen,qform,qwrite,ir)
         inquire(unit=ir,name=name,opened=qopen)
##IF PARALLEL
         CALL DDI_BCAST(1,'I',IR,1,0)
         CALL DDI_BCAST(1,'I',IW,1,0)
         call DDI_BCAST(1,'I',name,80,0)
##ENDIF
         write(*,*)'CHGMIU-a>ir,iw,q,n=',ir,iw,qopen,name(1:40)
      else
         IR=KSTRM
         IW=KOUTU
         CALL SEQOPN(IR,'INPUT', 'OLD',.TRUE., 'FORMATTED')
         CALL SEQOPN(IP,'PUNCH', 'NEW',.FALSE.,'FORMATTED')
         call getenv('INPUT',name)
         rewind(unit=ir)
         write(*,*)'CHGMIU-b>ir,iw,q,n=',ir,iw,qopen,name(1:40)
      endif
C
      RETURN
      END
C
      SUBROUTINE GAMEND(kstrm)
C-----------------------------------------------------------------------
C     Routine to put pointer in input file after all gamess input,
C     i.e. after last ' $end' string! 
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
CCC##INCLUDE '~/charmm_fcm/stream.fcm'
C
      integer kstrm
      CHARACTER*80 COMLYN
      INTEGER COMLEN,NENDS,IENDS
C
CCC      IF(IOLEV.GT.0) THEN
         REWIND KSTRM
         NENDS=0
 100     READ(KSTRM,45,ERR=945,END=955) COMLYN(1:80)
 45      FORMAT(A)
         COMLEN=80
         CALL CNVTUC(COMLYN,COMLEN)
         IF (INDX(COMLYN,COMLEN,'$END',4).GT.0) NENDS=NENDS+1
         GOTO 100
 955     CONTINUE
         WRITE(*,*)'NUMBER OF ENDs',NENDS
         REWIND KSTRM
         IENDS=0
 101     READ(KSTRM,45,ERR=945,END=945) COMLYN(1:80)
         COMLEN=80
         CALL CNVTUC(COMLYN,COMLEN)
         IF (INDX(COMLYN,COMLEN,'$END',4).GT.0) IENDS=IENDS+1
         IF(IENDS.EQ.NENDS) RETURN
         GOTO 101
         RETURN
 945     CALL WRNDIE(-2,'<GAMEND>','Unable to find $END')
CCC      ENDIF
C
      RETURN
      END
C
##IFN PARALLEL
C     The following is to enable compilation of CHARMM/GAMESS
C     without PARALLEL preflx keyword
C-----------------------------------------------------------------------
CCC      SUBROUTINE DDI_BCAST(MSGID,MSGTYP,BUFF,MSGLEN,IFROM)
CCC      DOUBLE PRECISION BUFF(*)
CCC      CHARACTER*1 MSGTYP
CCC      INTEGER MSGID,MSGLEN,IFROM
CCC      RETURN
CCC      END
      SUBROUTINE GCOMB(X,N)
      REAL*8 X(*)
      INTEGER N
      RETURN
      END
      INTEGER FUNCTION NNOD()
      NNOD=1
      RETURN
      END
      INTEGER FUNCTION MNOD()
      MNOD=0
      RETURN
      END
      SUBROUTINE PSYNC()
      RETURN
      END
##ENDIF
C
      SUBROUTINE COPSEL(ISLCT,QQINP)
C-----------------------------------------------------------------------
C     Copies selection vector to common block 
C     so it may be used by GAMESS interface
C     Call this routine only once and retain definition
C     of QM, MM, and link atoms throughout the claculation.
C     We call this from GAMINI which is called from charmm/charmm.src
C
C     IGMSEL(I) = 2  Link atom
C     IGMSEL(I) = 1  QM atom
C     IGMSEL(I) = 0  MM atom
C     IGMSEL(I) = -1 MM atom to be excluded from QM/MM interaction
C
C     MM atom in position close to link atom is excluded from interaction
C     of external charges to QM region. Instead of this atom is already
C     a link atom so no need for two atoms in one place!
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
C
      INTEGER ISLCT(*)
      LOGICAL QQINP
C
##INCLUDE '~/charmm_fcm/coord.fcm'
##INCLUDE '~/charmm_fcm/gamess.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/psf.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
C
      INTEGER I,J,I1,I2,N,LN,IS,IQ
      CHARACTER*4 SID, RID, REN, AC
      LOGICAL LNFLAG
C
      DO I=1, NATOM
         IGMSEL(I)=ISLCT(I)
         IF (TYPE(I)(1:2).EQ.'QQ') IGMSEL(I)=2
      ENDDO
C
C     Check if link atom is connected to any of its neighbors. If
C     yes then that atom will not be included in QM/MM interaction.
C     This is sometimes necessary to prevent oposite charge collision,
C     since QM cannot prevent this to happen.
C
C
      DO I=1,NBOND
         I1=IB(I)
         I2=JB(I)
         IF (IGMSEL(I1).EQ.2) THEN
C           Don't change QM atoms
            IF(QGMEXG) THEN
C              remove the entire group
               J=GETRES(I2,IGPBS,NGRP)
               IS=IGPBS(J)+1
               IQ=IGPBS(J+1)
CCC               write(*,*)'COPSEL>j,is,iq,i1,i2=',j,is,iq,i1,i2
               DO J=IS,IQ
                  IF(IGMSEL(J).EQ.0) IGMSEL(J)=-1
               ENDDO
            ELSE
C              remove the link host atom
               IF(IGMSEL(I2).EQ.0) IGMSEL(I2)=-1
            ENDIF
         ENDIF
         IF (IGMSEL(I2).EQ.2) THEN
            IF(QGMEXG) THEN
C              remove the entire group
               J=GETRES(I1,IGPBS,NGRP)
               IS=IGPBS(J)+1
               IQ=IGPBS(J+1)
CCC               write(*,*)'COPSEL>j,is,iq,i1,i2=',j,is,iq,i1,i2
               DO J=IS,IQ
                  IF(IGMSEL(J).EQ.0) IGMSEL(J)=-1
               ENDDO
            ELSE
C              remove the link host atom
               IF(IGMSEL(I1).EQ.0) IGMSEL(I1)=-1
            ENDIF
         ENDIF
      ENDDO
C
      IF(PRNLEV.GE.2) THEN
         WRITE(OUTU,118)
         WRITE(OUTU,120) 
     $        'Classical atoms excluded from the QM calculation'
         IF (QBLUCH) WRITE(OUTU,120) 
     $        'When BLUR option is used, it is controlled by WMAIN'
      ENDIF
 118  FORMAT('------------------------------------------------')
 120  FORMAT('GAMINT: ',A,':')
 122  FORMAT(10X,I5,4(1X,A4))
 124  FORMAT(10X,'NONE.')
      N=0
      DO I=1,NATOM
         IF(IGMSEL(I).EQ.-1) THEN
            CALL ATOMID(I,SID,RID,REN,AC)
            IF(PRNLEV.GE.2) WRITE(OUTU,122) I,SID,RID,REN,AC
            N=N+1
         ENDIF
      ENDDO
      IF(PRNLEV.GE.2) THEN
         IF(N.EQ.0) WRITE(OUTU,124)
         WRITE(OUTU,120) 'Quantum mechanical atoms'
      ENDIF
      N=0
      DO I=1,NATOM
         IF(IGMSEL(I).EQ.1) THEN
            CALL ATOMID(I,SID,RID,REN,AC)
            IF(PRNLEV.GE.2) WRITE(OUTU,122) I,SID,RID,REN,AC
            N=N+1
         ENDIF
      ENDDO
      IF(PRNLEV.GE.2) THEN
         IF(N.EQ.0) WRITE(OUTU,124)
         WRITE(OUTU,120) 'Quantum mechanical link atoms'
      ENDIF
      N=0
      DO I=1,NATOM
         IF(IGMSEL(I).EQ.2) THEN
            CALL ATOMID(I,SID,RID,REN,AC)
            IF(PRNLEV.GE.2) WRITE(OUTU,122) I,SID,RID,REN,AC
            N=N+1
         ENDIF
      ENDDO
      IF(PRNLEV.GE.2) THEN
         IF(N.EQ.0) WRITE(OUTU,124)
         WRITE(OUTU,118)
      ENDIF
C
C     Alow for partial charges on any QM atom
C
      N=0
      LN=0
      DO I = 1,NATOM
        IF (IGMSEL(I) .GE. 1) THEN
           N = N + 1
C
C     Non integer charges for link atoms can be specified separately
C     Also allow to change them subsequently with SCALar command 
C     using QINP keyword in GAMEss command. 
C
           LNFLAG=.FALSE.
           IF((TYPE(I)(1:2).EQ.'QQ').AND.(NQQCHG.NE.0)) THEN
              LN=LN+1
              IF(QQCHG(LN).GT.-NINE99) THEN
                 FQQCHG(N) = QQCHG(LN)
              ELSE
                 FQQCHG(N) = -THOSND
              ENDIF
C
C     Don't have any more link atoms, put NQQCHG to 0
C     for possible subsequent changes with SCALar command
C     or on restarts.
C
              IF(LN.EQ.NQQCHG) NQQCHG=0
              LNFLAG=.TRUE.
           ENDIF
C
C
C     First check the flag QINP for all atoms.
C
           IF(QQINP) THEN
C
C     All QM charges (accept ones specified in ADDLink
C     are taken from PSF
C
              IF(.NOT.LNFLAG) FQQCHG(N)=CG(I)
           ELSE
              FQQCHG(N)=-THOSND
           ENDIF
C
           IF(PRNLEV.GE.2) WRITE(OUTU,'(A,2I5,A,F15.5)')
     $          'GAMINT: ATOM(',I,N,') has QNUC: ',FQQCHG(N)
        ENDIF
      ENDDO
C
C Zero charges on quantum atoms to remove from MM term.
      IF(QGMREM) THEN
         DO I=1, NATOM
            IF(IGMSEL(I).GT.0) CG(I)=ZERO
         ENDDO
      ENDIF
CCC      write(*,*)'COPSEL>IGMSEL=',(IGMSEL(I),I=1,NATOM)
C
##ENDIF
##ENDIF
      RETURN
      END

##IF GAMESS
##IFN FLUCQ

      SUBROUTINE FQQMMM(DMAT)
C Dummy subroutine to ensure GAMESS-US links properly with CHARMM when
C FlucQ is not being used
      DOUBLE PRECISION DMAT(*)
      RETURN
      END

##ENDIF
##ENDIF
