text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   Potential energy routine taken from the POTLIB library:
!   R. J. Duchovic, Y. L. Volobuev, G. C. Lynch, A. W. Jasper, D. G. Truhlar, T. C. Allison,
!   A. F. Wagner, B. C. Garrett, J. Espinosa-García, and J. C. Corchado, POTLIB,
!   http://comp.chem.umn.edu/potlib.

c
c     driver routine
c
      subroutine egrad_clnh3(q,Natoms,Nbeads,V,dVdq,info)
      integer, intent(in) :: Natoms
      integer, intent(in) :: Nbeads
      double precision, intent(in) :: q(3,Natoms,Nbeads)
      double precision, intent(out) :: V(Nbeads)
      double precision, intent(out) :: dVdq(3,Natoms,Nbeads)
      double precision :: R(15), DEGSDR(15),ENGYGS
      integer::i,j

c
c     convert input structure
c
      do i=1,natoms
         do j=1,3
            R((i-1)*3+j)=q(j,i,1)
         end do
      end do
      call pot_clnh3(R,ENGYGS,DEGSDR)

c
c     convert output gradient
c
      V(1)=ENGYGS
      do i=1,natoms
         do j=1,3
            dVdq(j,i,1)=DEGSDR((i-1)*3+j)
         end do
      end do
c      stop "Jjgfdf"
      return
      end subroutine

c**************************************************************************
C   System:           ClNH3
C   Number of electronic surfaces: 1
C   Number of derivatives: 0
C   Number of bodies: 5
C   Interface:        potlib2001
C   Common name:      ClNH3
C   Notes:
C
C   References: M. Monge-Palacios, C. Rangel, J.C. Corchado and J. Espinosa-Garcia,
C   Int.J.Quantum.Chem. 112, 1887 (2012)

      SUBROUTINE POT_clnh3(R_out,ENGYGS_out,DEGSDR_out)
C
C        This potential is written such that:
C
C                       X(1)  - X(3)  : X, Y, Z for H1
C                       X(4)  - X(6)  : X, Y, Z for N
C                       X(7)  - X(9)  : X, Y, Z for H2
C                       X(10) - X(12) : X, Y, Z for H3
C                       X(13) - X(15) : X, Y, Z for Cl (Hb)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      real(kind=8)::R_out(15),ENGYGS_out,DEGSDR_out(15)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
C      DIMENSION COORD(N3TM),DX(N3TM)
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
      CALL CARTOU_clnh3
      CALL CARTTOR_clnh3
C
C Changing units from angstrom to bohr and initialize
C
      R(1:15)=R_OUT
      DO 10 I = 1, 15
         q(I) = R(I)*0.52918d0
         pdot(I)=0.D0
 10   CONTINUE
c
c  calculate relative coordinates and bond lengths
c
       en=0.0d0

       call coorden_clnh3
c
c  calculate switching functions
c
       call switchf_clnh3
c
c  calculate reference angles and their derivatives
c
       call refangles_clnh3
       call oprefangles_clnh3
c
c  calculate stretching potential
c
       call stretch_clnh3(vstr)
c
c  calculate out of plane bending potential
c
c      call opbend(vop)
c
c  calculate in plane bending potential
c
       call ipbend_clnh3(vip)
c
c  total potential energy is vstr+vop+vip
c
       en=vstr+vip
c      en=vstr+vop+vip
c
c 0.03812 conversion factor from 10(5) j/mol to hartrees
c
       en = en*0.03812D0
       ENGYGS = en
c
      CALL EUNITZERO_clnh3
      IF(NDER.NE.0) THEN
c
c 0.0201723 conversion factor from 10(5)j/mol/A to hartrees/bohr
c
         do i=1,15
            DEGSDR(i)=pdot(i)*0.0201723d0
         enddo
         CALL RTOCART_clnh3
         CALL DEDCOU_clnh3
      ENDIF
      ENGYGS_OUT=ENGYGS
      DEGSDR_OUT=DEGSDR(1:15)

C
       return
       end
c
c******************************************************
c
       subroutine coorden_clnh3
c
c  calculates relative coordinates and bond lengths
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
c  calculate relative coordinates
c
       do ind=1,3
         tcb(ind)=q(nc(ind))-q(nhb(ind))
         do i=1,3
           tch(i,ind)=q(nc(ind))-q(nh(i,ind))
           tbh(i,ind)=q(nhb(ind))-q(nh(i,ind))
         enddo
       enddo
c
c  calculate bond lengths
c
       rcb=sqrt(tcb(1)*tcb(1)+tcb(2)*tcb(2)+tcb(3)*tcb(3))
       do i=1,3
         rch(i)=sqrt(tch(i,1)*tch(i,1)+tch(i,2)*tch(i,2)+
     *                tch(i,3)*tch(i,3))
         rbh(i)=sqrt(tbh(i,1)*tbh(i,1)+tbh(i,2)*tbh(i,2)+
     *                tbh(i,3)*tbh(i,3))
       enddo
c jcc-2010
       argmax = 19.d0
       P1 = 1.d0
       do i = 1, 3
          argp1 = (w1*(rch(i)-w2))
          if (argp1.lt.argmax) then
            t1tmp = 1 - tanh(argp1)
          else
            t1tmp = 0.d0
          endif
          P1 = P1 * t1tmp
       enddo
       r0ch = P1 * r0chr + (1-P1) * r0chp
c jcc-2010
       return
       end
c
c******************************************************
c
c
       subroutine refangles_clnh3
c
c  subroutine calculates reference angles for the "in-plane" potential
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
C      tau=acos(-1.0d0/3.0d0)
C      pi=4.0d0*atan(1.0d0)
       halfpi=0.5d0*pi
       twopi=2.0d0*pi
       ppito=(twopi-taunh2)/2.0d0
c
c  set diagonal elements to zero
c
       do i=1,3
         theta0(i,i)=0.0d0
         do k=1,4
           dtheta0(i,i,k)=0.0d0
         enddo
       enddo
c
c  calculate reference angles
c
       theta0(1,2)=tau+(tau-ppito)*(sphi(1)*sphi(2)-1.0d0)
     *             +(tau-taunh2)*(stheta(3)-1.0d0)
       theta0(1,3)=tau+(tau-ppito)*(sphi(1)*sphi(3)-1.0d0)
     *             +(tau-taunh2)*(stheta(2)-1.0d0)
       theta0(2,3)=tau+(tau-ppito)*(sphi(2)*sphi(3)-1.0d0)
     *             +(tau-taunh2)*(stheta(1)-1.0d0)
c
c  calculate the derivatives of theta0(i,j) in terms of rch(k)
c  quantity calulated is dtheta0(i,j,k)
c
c  derivatives wrt rch(1)
c
       dtheta0(1,2,1)=(tau-ppito)*dsphi(1)*sphi(2)
       dtheta0(1,3,1)=(tau-ppito)*dsphi(1)*sphi(3)
       dtheta0(2,3,1)=(tau-taunh2)*dstheta(1)
c
c  derivatives wrt rch(2)
c
       dtheta0(1,2,2)=(tau-ppito)*sphi(1)*dsphi(2)
       dtheta0(1,3,2)=(tau-taunh2)*dstheta(2)
       dtheta0(2,3,2)=(tau-ppito)*dsphi(2)*sphi(3)
c
c  derivatives wrt rch(3)
c
       dtheta0(1,2,3)=(tau-taunh2)*dstheta(3)
       dtheta0(1,3,3)=(tau-ppito)*sphi(1)*dsphi(3)
       dtheta0(2,3,3)=(tau-ppito)*sphi(2)*dsphi(3)
c
c  fill in the other half of the matrix
c
        do i=1,2
          do j=i+1,3
            theta0(j,i)=theta0(i,j)
            do k=1,3
              dtheta0(j,i,k)=dtheta0(i,j,k)
            enddo
          enddo
        enddo
       return
       end
c
c******************************************************
c
c
       subroutine oprefangles_clnh3
c
c  subroutine calculates reference angles for the "out-plane" potential
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
C       pi=4.0d0*atan(1.0d0)
       halfpi=0.5d0*pi
       twopi=2.0d0*pi
       ppito=(twopi-taunh2)/2.0d0
c
c  set diagonal elements to zero
c
       do i=1,2
         optheta0(i,i)=0.0d0
         do k=1,3
           doptheta0(i,i,k)=0.0d0
         enddo
       enddo
c
c  calculate reference angles
c
       optheta0(1,1)=optau+(optau-halfpi)*(sphi(1)
     *          -1.0d0)+(optau-halfpi)*(stheta(3)-1.0d0)
     *          +(optau-halfpi)*(stheta(2)-1.0d0)
       optheta0(2,2)=optau+(optau-halfpi)*(sphi(2)
     *          -1.0d0)+(optau-halfpi)*(stheta(3)-1.0d0)
     *          +(optau-halfpi)*(stheta(1)-1.0d0)
       optheta0(3,3)=optau+(optau-halfpi)*(sphi(3)
     *          -1.0d0)+(optau-halfpi)*(stheta(1)-1.0d0)
     *          +(optau-halfpi)*(stheta(2)-1.0d0)
c
       optheta0(1,2)=optau+(optau-halfpi)*(sphi(1)*sphi(2)
     *          -1.0d0)+(optau-halfpi)*(stheta(3)-1.0d0)
       optheta0(1,3)=optau+(optau-halfpi)*(sphi(1)*sphi(3)
     *          -1.0d0)+(optau-halfpi)*(stheta(2)-1.0d0)
       optheta0(2,3)=optau+(optau-halfpi)*(sphi(2)*sphi(3)
     *          -1.0d0)+(optau-halfpi)*(stheta(1)-1.0d0)
c
c  calculate the derivatives of optheta0(i,j) in terms of rch(k)
c  quantity calulated is doptheta0(i,j,k)
c
c  derivatives wrt rch(1)
c
       doptheta0(1,1,1)=(optau-halfpi)*dsphi(1)
       doptheta0(2,2,1)=(optau-halfpi)*dstheta(1)
       doptheta0(3,3,1)=(optau-halfpi)*dstheta(1)
c
       doptheta0(1,2,1)=(optau-halfpi)*dsphi(1)*
     * sphi(2)
       doptheta0(1,3,1)=(optau-halfpi)*dsphi(1)*
     * sphi(3)
       doptheta0(2,3,1)=(optau-halfpi)*
     * dstheta(1)
c
c  derivatives wrt rch(2)
c
       doptheta0(1,1,2)=(optau-halfpi)*dstheta(2)
       doptheta0(2,2,2)=(optau-halfpi)*dsphi(2)
       doptheta0(3,3,2)=(optau-halfpi)*dstheta(2)
c
       doptheta0(1,2,2)=(optau-halfpi)*sphi(1)*
     * dsphi(2)
       doptheta0(1,3,2)=(optau-halfpi)*
     * dstheta(2)
       doptheta0(2,3,2)=(optau-halfpi)*dsphi(2)*
     * sphi(3)
c
c  derivatives wrt rch(3)
c
c
       doptheta0(1,1,3)=(optau-halfpi)*dstheta(3)
       doptheta0(2,2,3)=(optau-halfpi)*dstheta(3)
       doptheta0(3,3,3)=(optau-halfpi)*dsphi(3)
c
       doptheta0(1,2,3)=(optau-halfpi)*
     * dstheta(3)
       doptheta0(1,3,3)=(optau-halfpi)*sphi(1)*
     * dsphi(3)
       doptheta0(2,3,3)=(optau-halfpi)*sphi(2)*
     * dsphi(3)
c
c
c  fill in the other half of the matrix
c
        do i=1,2
          do j=i+1,3
            optheta0(j,i)=optheta0(i,j)
            do k=1,3
              doptheta0(j,i,k)=doptheta0(i,j,k)
            enddo
          enddo
        enddo
       return
       end
c
c******************************************************
c
c
       subroutine stretch_clnh3(vstr)
c
c  subroutine to calculate leps-type stretching potential and its
c  derivatives
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       dimension vqch(4),vjch(4),vqbh(4),vjbh(4),vq(4),vj(4),
     *           achdc(3),achdh(4,3)
c
c  calculate avergage bond length for the nh3 moiety
c
       rav=(rch(1)+rch(2)+rch(3))/3.0d0
c
c  initialise:
c
       vstr=0.0d0
c
c  ach:
c
c  in double precision tanh(19.0d0)=1.0d0 and we put the if statement
c  in to avoid overflow/underflow errors
c
       arga=c1ch*(rav-r0ch)
       if(arga.lt.19.0d0)then
         ach=a1ch+b1ch*(tanh(arga)+1.0d0)*0.5d0
         dumach=b1ch*c1ch/(2.0d0*cosh(arga)**2)
       else
         ach=a1ch+b1ch
         dumach=0.0d0
       endif
c
c  calculate singlet: e1, triplet: e3 energies and vq and vj
c  terms for each bond
c
       e1=d1cb*(exp(-2.0d0*acb*(rcb-r0cb))-2.0d0*exp(-acb*(rcb-r0cb)))
       e3=d3cb*(exp(-2.0d0*acb*(rcb-r0cb))+2.0d0*exp(-acb*(rcb-r0cb)))
       vqcb=(e1+e3)*0.5d0
       vjcb=(e1-e3)*0.5d0
       do i=1,3
         e1=d1ch*(exp(-2.0d0*ach*(rch(i)-r0ch))
     *              -2.0d0*exp(-ach*(rch(i)-r0ch)))
         e3=d3ch*(exp(-2.0d0*ach*(rch(i)-r0ch))
     *              +2.0d0*exp(-ach*(rch(i)-r0ch)))
         vqch(i)=(e1+e3)*0.5d0
         vjch(i)=(e1-e3)*0.5d0
         e1=d1hh*(exp(-2.0d0*ahh*(rbh(i)-r0hh))
     *              -2.0d0*exp(-ahh*(rbh(i)-r0hh)))
         e3=d3hh*(exp(-2.0d0*ahh*(rbh(i)-r0hh))
     *              +2.0d0*exp(-ahh*(rbh(i)-r0hh)))
         vqbh(i)=(e1+e3)*0.5d0
         vjbh(i)=(e1-e3)*0.5d0
c
c  calculate 3 body potential
c
         vq(i)=vqch(i)+vqcb+vqbh(i)
         vj(i)=-sqrt(((vjch(i)-vjcb)**2+(vjcb-vjbh(i))**2
     *                 +(vjbh(i)-vjch(i))**2)*0.5d0)
         vstr=vstr+vq(i)+vj(i)
       enddo
c
c  partial derivatives
c  first we need the derivative of ach:
c
       do ind=1,3
         achdc(ind)=dumach*(tch(1,ind)/rch(1)+tch(2,ind)/rch(2)
     *            +tch(3,ind)/rch(3))/3.0d0
         do i=1,3
           achdh(i,ind)=-dumach*tch(i,ind)/rch(i)/3.0d0
         enddo
       enddo
       dumqcb=-acb*((d1cb+d3cb)*exp(-2.0d0*acb*(rcb-r0cb))-
     *         (d1cb-d3cb)*exp(-acb*(rcb-r0cb)))/rcb
c
c  calculate cartesian derivatives:
c  looping over ch(i) and bh(i)
c
       do i=1,3
         dumqbh=-ahh*((d1hh+d3hh)*exp(-2.0d0*ahh*(rbh(i)-r0hh))-
     *           (d1hh-d3hh)*exp(-ahh*(rbh(i)-r0hh)))/rbh(i)
         factj=0.5d0/vj(i)
         dumjcb=-acb*((d1cb-d3cb)*exp(-2.0d0*acb*(rcb-r0cb))
     *            -(d1cb+d3cb)*exp(-acb*(rcb-r0cb)))*factj/rcb
         dumjbh=-ahh*((d1hh-d3hh)*exp(-2.0d0*ahh*(rbh(i)-r0hh))
     *            -(d1hh+d3hh)*exp(-ahh*(rbh(i)-r0hh)))*factj/rbh(i)
         do ind=1,3
c
c  deriv wrt hb:
c
                  pdot(nhb(ind))=pdot(nhb(ind))
     *             -tcb(ind)*dumqcb+tbh(i,ind)*dumqbh
     *            +(vjch(i)-vjcb)*(dumjcb*tcb(ind))
     *            +(vjcb-vjbh(i))*(-dumjcb*tcb(ind)-dumjbh*tbh(i,ind))
     *            +(vjbh(i)-vjch(i))*dumjbh*tbh(i,ind)
c
c  dvqch(i)/dc
c
           dumqch=-(ach*tch(i,ind)/rch(i)+achdc(ind)*(rch(i)-r0ch))
     *              *((d1ch+d3ch)*exp(-2.0d0*ach*(rch(i)-r0ch))
     *                 -(d1ch-d3ch)*exp(-ach*(rch(i)-r0ch)))
               pdot(nc(ind))=pdot(nc(ind))+dumqch+tcb(ind)*dumqcb
c
c  dvqch(i)/dh(i)
c
           dumqhi=(ach*tch(i,ind)/rch(i)-achdh(i,ind)*(rch(i)-r0ch))
     *              *((d1ch+d3ch)*exp(-2.0d0*ach*(rch(i)-r0ch))
     *                 -(d1ch-d3ch)*exp(-ach*(rch(i)-r0ch)))
              pdot(nh(i,ind))=pdot(nh(i,ind))+dumqhi-tbh(i,ind)*dumqbh
c
c  dvjch(i)/dc
c
           dumjch=-(ach*tch(i,ind)/rch(i)+achdc(ind)*(rch(i)-r0ch))
     *              *((d1ch-d3ch)*exp(-2.0d0*ach*(rch(i)-r0ch))
     *               -(d1ch+d3ch)*exp(-ach*(rch(i)-r0ch)))*factj
c
c  dvj(i)/dnc(ind)
c
           pdot(nc(ind))=pdot(nc(ind))
     *            +(vjch(i)-vjcb)*(dumjch-dumjcb*tcb(ind))
     *            +(vjcb-vjbh(i))*dumjcb*tcb(ind)
     *            -(vjbh(i)-vjch(i))*dumjch
c
c  dvjch(i)/dh(i)
c
           dumjhi=(ach*tch(i,ind)/rch(i)-achdh(i,ind)*(rch(i)-r0ch))
     *              *((d1ch-d3ch)*exp(-2.0d0*ach*(rch(i)-r0ch))
     *               -(d1ch+d3ch)*exp(-ach*(rch(i)-r0ch)))*factj
c
c  dvj(i)/dnh(i,ind)
c
            pdot(nh(i,ind))=pdot(nh(i,ind))
     *            +(vjch(i)-vjcb)*dumjhi
     *            +(vjcb-vjbh(i))*dumjbh*tbh(i,ind)
     *            +(vjbh(i)-vjch(i))*(-dumjbh*tbh(i,ind)-dumjhi)
c
c  dv(i)/dh(j)
c
           do k=1,2
             j=i+k
             if(j.gt.3)j=j-3
             dumqhj=-achdh(j,ind)*(rch(i)-r0ch)
     *                 *((d1ch+d3ch)*exp(-2.0d0*ach*(rch(i)-r0ch))
     *                    -(d1ch-d3ch)*exp(-ach*(rch(i)-r0ch)))
             dumjhj=-achdh(j,ind)*(rch(i)-r0ch)
     *                 *((d1ch-d3ch)*exp(-2.0d0*ach*(rch(i)-r0ch))
     *                  -(d1ch+d3ch)*exp(-ach*(rch(i)-r0ch)))*factj
             pdot(nh(j,ind))=pdot(nh(j,ind))+dumqhj
     *            +(vjch(i)-vjcb)*dumjhj
     *            -(vjbh(i)-vjch(i))*dumjhj
           enddo
         enddo
       enddo
       return
       end
c
c******************************************************
c
c
       subroutine opbend_clnh3(vop)
c
c  subroutine calculates symmetrized vop potential and derivatives
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       double precision norma
       dimension sumd2(4),sumd4(4)
       dimension in(3),a(3),b(3),axb(3),c(4,3),argd(4)
c
c
       vop=0.0d0
c
c  calculate force constants and their derivatives
c
       call opforce
c
c  calculate out-of-plane angle and derivatives
c
       do i=1,3
         j=i+1
         if(j.gt.3)j=j-3
         k=j+1
         if(k.gt.3)k=k-3
c
c  modification to ensure that the set of methane CH bond vector
c  (rj,rk,rl) is a right-handed set
c
       in(1)=i
       in(2)=j
       in(3)=k
c
c  vector a is rk-rj, vector b is rl-rj
c
       do ind=1,3
         a(ind)=q(nh(j,ind))-q(nh(i,ind))
         b(ind)=q(nh(k,ind))-q(nh(i,ind))
       enddo
c
c  axb is vector a cross b
c
       axb(1)=a(2)*b(3)-a(3)*b(2)
       axb(2)=a(3)*b(1)-a(1)*b(3)
       axb(3)=a(1)*b(2)-a(2)*b(1)
       norma=axb(1)*axb(1)+axb(2)*axb(2)+axb(3)*axb(3)
       norma=sqrt(norma)
c
c  c is position vector of h(ii): calculate c(j),c(k),c(l)
c
       do ii=1,3
         do ind=1,3
           c(in(ii),ind)=-tch(in(ii),ind)/rch(in(ii))
         enddo
       enddo
c
c  argd is the dot product axb dot c
c
       do ii=1,3
         argd(in(ii))=axb(1)*c(in(ii),1)+axb(2)*c(in(ii),2)
     *                                +axb(3)*c(in(ii),3)
         argd(in(ii))=argd(in(ii))/norma
c
c  if argd < 0 we need to switch vectors k and l around
c
         if (argd(in(ii)).lt.0.d0) then
             itemp=j
             j=k
             k=itemp
         endif
       enddo
c
c  sum2 = sum delta**2
c  sum4 = sum delta**4
c
         call calcdelta(i,j,k,sum2,sum4)
         sumd2(i)=sum2
         sumd4(i)=sum4
         vop=vop+fdelta(i)*sumd2(i)+hdelta(i)*sumd4(i)
       enddo
c
       do i=1,3
         do j=1,3
c
c  overall derivatives of force constants i wrt the bond-length rch(j)
c
           ddr=dfdelta(i,j)*sumd2(i)+dhdelta(i,j)*sumd4(i)
c
c  calculate derivatives in terms of cartesian coordinates:
c
           do ind=1,3
             pdot(nh(j,ind))=pdot(nh(j,ind))-tch(j,ind)*ddr/rch(j)
             pdot(nc(ind))=pdot(nc(ind))+tch(j,ind)*ddr/rch(j)
           enddo
         enddo
       enddo
       return
       end
c
c******************************************************
c
c
       subroutine ipbend_clnh3(vip)
c
c  subroutine calculates symmetrised in plane bend term
c  and its derivatives
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       dimension costh(4,4),theta(4,4),dth(4,4)
c
c  initialise
c
       vip=0.0d0
c
c  calculate force constants: fk0(i,j), f1(i)
c  and derivatives wrt rch(k) and rbh(k): dfdc(i,j,k), dfdh(i,j,k)
c
       call ipforce_clnh3
c
c  calculate theta(i,j) and in plane bend potential
c
       do i=1,2
         do j=i+1,3
           costh(i,j)=tch(i,1)*tch(j,1)+tch(i,2)*tch(j,2)
     *                       +tch(i,3)*tch(j,3)
           costh(i,j)=costh(i,j)/rch(i)/rch(j)
           theta(i,j)=acos(costh(i,j))
           dth(i,j)=theta(i,j)-theta0(i,j)
           vip=vip+0.5d0*fk0(i,j)*f1(i)*f1(j)*dth(i,j)**2
c
c  calculate partial derivatives wrt cartesian coordinates
c
c  calculate pdots wrt theta:
c
           termth=-1.0d0/sqrt(1.0d0-costh(i,j)*costh(i,j))
           do ind=1,3
             dthi=-tch(j,ind)/rch(i)/rch(j)
     *                  +costh(i,j)*tch(i,ind)/rch(i)/rch(i)
             dthi=dthi*termth
             dthj=-tch(i,ind)/rch(i)/rch(j)
     *                  +costh(i,j)*tch(j,ind)/rch(j)/rch(j)
             dthj=dthj*termth
             dthc=-(dthi+dthj)
             pdot(nh(i,ind))=pdot(nh(i,ind))
     *                       +fk0(i,j)*f1(i)*f1(j)*dthi*dth(i,j)
             pdot(nh(j,ind))=pdot(nh(j,ind))
     *                       +fk0(i,j)*f1(i)*f1(j)*dthj*dth(i,j)
             pdot(nc(ind))=pdot(nc(ind))
     *                       +fk0(i,j)*f1(i)*f1(j)*dthc*dth(i,j)
             do k=1,3
c
c  calculate pdots wrt force constants and wrt theta0
c
               dth0k=-dtheta0(i,j,k)*tch(k,ind)/rch(k)
               dth0c=-dth0k
               pdot(nh(k,ind))=pdot(nh(k,ind))
     *                -0.5d0*tch(k,ind)*dfdc(i,j,k)*dth(i,j)**2/rch(k)
     *                -0.5d0*tbh(k,ind)*dfdh(i,j,k)*dth(i,j)**2/rbh(k)
     *                      -fk0(i,j)*f1(i)*f1(j)*dth0k*dth(i,j)
               pdot(nc(ind))=pdot(nc(ind))
     *                +0.5d0*tch(k,ind)*dfdc(i,j,k)*dth(i,j)**2/rch(k)
     *                      -fk0(i,j)*f1(i)*f1(j)*dth0c*dth(i,j)
               pdot(nhb(ind))=pdot(nhb(ind))
     *                +0.5d0*tbh(k,ind)*dfdh(i,j,k)*dth(i,j)**2/rbh(k)
             enddo
           enddo
         enddo
       enddo
       return
       end
c
c*************************************************************************
c
       subroutine calcdelta_clnh3(i,j,k,sum2,sum4)
c
c  subroutine calculates out of plane angle delta, loops
c  through delta(i,j), delta(i,k), delta(i,l)
c
c   also calculates the derivatives wrt delta
c
       implicit double precision (a-h,o-z)
       double precision norma
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       dimension  delta(4),in(3),a(3),b(3),axb(3),c(4,3),argd(4),
     *            daxb(4,3,3),cdot(4,3,3),atemp2(3)
c
c  set i,j,k indices
c
       in(1)=i
       in(2)=j
       in(3)=k
c
c  vector a is rj-ri, vector b is rk-ri
c
       do ind=1,3
         a(ind)=q(nh(j,ind))-q(nh(i,ind))
         b(ind)=q(nh(k,ind))-q(nh(i,ind))
       enddo
c
c  axb is vector a cross b
c
       axb(1)=a(2)*b(3)-a(3)*b(2)
       axb(2)=a(3)*b(1)-a(1)*b(3)
       axb(3)=a(1)*b(2)-a(2)*b(1)
       norma=axb(1)*axb(1)+axb(2)*axb(2)+axb(3)*axb(3)
       norma=sqrt(norma)
c
c  c is position vector of h(ii): calculate c(j),c(k),c(l)
c
       do ii=1,3
         do ind=1,3
           c(in(ii),ind)=-tch(in(ii),ind)/rch(in(ii))
         enddo
       enddo
c
c  initialise
c
       sum2=0.0d0
       sum4=0.0d0
c
c  argd is the dot product axb dot c
c
       do ii=1,3
         argd(in(ii))=axb(1)*c(in(ii),1)+axb(2)*c(in(ii),2)
     *                                +axb(3)*c(in(ii),3)
         argd(in(ii))=argd(in(ii))/norma
         delta(in(ii))=acos(argd(in(ii)))-optheta0(i,in(ii))
c
         sum2=sum2+delta(in(ii))**2
         sum4=sum4+delta(in(ii))**4
       enddo
c
c  derivatives of axb wrt hi:
c
       daxb(i,1,1)=0.0d0
       daxb(i,1,2)=b(3)-a(3)
       daxb(i,1,3)=-b(2)+a(2)
       daxb(i,2,1)=-b(3)+a(3)
       daxb(i,2,2)=0.0d0
       daxb(i,2,3)=b(1)-a(1)
       daxb(i,3,1)=b(2)-a(2)
       daxb(i,3,2)=-b(1)+a(1)
       daxb(i,3,3)=0.0d0
c
c  derivatives of axb wrt hj:
c
       daxb(j,1,1)=0.0d0
       daxb(j,1,2)=-b(3)
       daxb(j,1,3)=b(2)
       daxb(j,2,1)=b(3)
       daxb(j,2,2)=0.0d0
       daxb(j,2,3)=-b(1)
       daxb(j,3,1)=-b(2)
       daxb(j,3,2)=b(1)
       daxb(j,3,3)=0.0d0
c
c  derivatives of axb wrt hl:
c
       daxb(k,1,1)=0.0d0
       daxb(k,1,2)=a(3)
       daxb(k,1,3)=-a(2)
       daxb(k,2,1)=-a(3)
       daxb(k,2,2)=0.0d0
       daxb(k,2,3)=a(1)
       daxb(k,3,1)=a(2)
       daxb(k,3,2)=-a(1)
       daxb(k,3,3)=0.0d0
c
c   loop over cdot(in(ii),ind,jind) where we consider deriv of c(in(ii))
c   wrt h(in(ii),jind) with components jind
c
       do ii=1,3
c
c  deriv of cdot(in(ii),x) wrt x, y, z
c
         cdot(in(ii),1,1)=1.0d0/rch(in(ii))
     *                   +tch(in(ii),1)*c(in(ii),1)/rch(in(ii))**2
         cdot(in(ii),1,2)=tch(in(ii),2)*c(in(ii),1)/rch(in(ii))**2
         cdot(in(ii),1,3)=tch(in(ii),3)*c(in(ii),1)/rch(in(ii))**2
c
c  deriv of cdot(in(ii),y) wrt x, y, z
c
         cdot(in(ii),2,1)=tch(in(ii),1)*c(in(ii),2)/rch(in(ii))**2
         cdot(in(ii),2,2)=1.0d0/rch(in(ii))
     *                   +tch(in(ii),2)*c(in(ii),2)/rch(in(ii))**2
         cdot(in(ii),2,3)=tch(in(ii),3)*c(in(ii),2)/rch(in(ii))**2
c
c  deriv of cdot(in(ii),z) wrt x, y, z
c
         cdot(in(ii),3,1)=tch(in(ii),1)*c(in(ii),3)/rch(in(ii))**2
         cdot(in(ii),3,2)=tch(in(ii),2)*c(in(ii),3)/rch(in(ii))**2
         cdot(in(ii),3,3)=1.0d0/rch(in(ii))
     *                   +tch(in(ii),3)*c(in(ii),3)/rch(in(ii))**2
       enddo
c
       do ii=1,3
         do ind=1,3
            deldot=-doptheta0(i,ii,i)
c
c  derivative wrt h(i,ind)
c  for  rch(i) only terms are from the derivatives of theta0
c
            deldot=-deldot*tch(i,ind)/rch(i)
c
c  derivative wrt c(ind)
c
            deldot=-deldot
c
           do jj=1,3
c
c  partial derivatives wrt h(in(jj),ind), loop over delta(i,in(ii))
c
c   atemp1 is axb dot daxb wrt h(in(jj))
c
            atemp1=axb(1)*daxb(in(jj),ind,1)
     *            +axb(2)*daxb(in(jj),ind,2)
     *            +axb(3)*daxb(in(jj),ind,3)
            atemp1=atemp1/(norma**3)
c
c  atemp2 is deriv of normalised axb
c
            atemp2(1)=daxb(in(jj),ind,1)/norma-atemp1*axb(1)
            atemp2(2)=daxb(in(jj),ind,2)/norma-atemp1*axb(2)
            atemp2(3)=daxb(in(jj),ind,3)/norma-atemp1*axb(3)
c
c  atemp3 is daxb dot c(in(ii))

            atemp3=atemp2(1)*c(in(ii),1)+atemp2(2)*c(in(ii),2)
     *                             +atemp2(3)*c(in(ii),3)
c
c  atemp4 is axb dot cdot
c
            atemp4=0.0d0
            if(ii.eq.jj)then
c
c  ie deriv of c(in(ii)) wrt h(in(jj)) is non zero only for ii = jj
c
              atemp4=axb(1)*cdot(in(ii),1,ind)
     *                     +axb(2)*cdot(in(ii),2,ind)
     *                     +axb(3)*cdot(in(ii),3,ind)
              atemp4=atemp4/norma
            endif
c
c  atemp5 is deriv of theta0(i,in(ii)) wrt to nh(in(jj),ind)
c
            atemp5=-doptheta0(i,in(ii),in(jj))
c
c  deriv wrt h(in(jj)),ind):
c
            atemp5=-atemp5*tch(in(jj),ind)/rch(in(jj))
            deldot=atemp3+atemp4
            deldot=-1.0d0/sqrt(1.0d0-argd(in(ii))**2)*deldot
            deldot=deldot+atemp5
            pdot(nh(in(jj),ind))=pdot(nh(in(jj),ind))
     *                   +2.0d0*fdelta(i)*delta(in(ii))*deldot
     *                   +4.0d0*hdelta(i)*delta(in(ii))**3*deldot
c
c  for carbon the only contributions are from axb dot cdot term and
c  from theta0 and derivative cdot wrt carbon=-cdot wrt hydrogen
c
            deldot=1.0d0/sqrt(1.0d0-argd(in(ii))**2)*atemp4
            deldot=deldot-atemp5
            pdot(nc(ind))=pdot(nc(ind))
     *                   +2.0d0*fdelta(i)*delta(in(ii))*deldot
     *                   +4.0d0*hdelta(i)*delta(in(ii))**3*deldot
          enddo
        enddo
       enddo
       return
       end

c******************************************************
c
       subroutine opforce_clnh3
c
c  calculates the out-of-plane bending force constants
c  and their derivatives
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       dimension switch(4),dswitch(4,4)
c
c  calculate switching functions:
c
       switch(1)=s3(1)*s3(2)*s3(3)
       switch(2)=s3(2)*s3(3)*s3(1)
       switch(3)=s3(3)*s3(1)*s3(2)
c
c  calculate derivatives:
c  derivative of switch(1) wrt the 3 rch bond lengths
c
       dswitch(1,1)=ds3(1)*s3(2)*s3(3)
       dswitch(1,2)=s3(1)*ds3(2)*s3(3)
       dswitch(1,3)=s3(1)*s3(2)*ds3(3)
c
c  derivative of switch(2) wrt the 3 rch bond lengths
c
       dswitch(2,1)=s3(2)*s3(3)*ds3(1)
       dswitch(2,2)=ds3(2)*s3(3)*s3(1)
       dswitch(2,3)=s3(2)*ds3(3)*s3(1)
c
c  derivative of switch(3) wrt the 3 rch bond lengths
c
       dswitch(3,1)=s3(3)*ds3(1)*s3(2)
       dswitch(3,2)=s3(3)*s3(1)*ds3(2)
       dswitch(3,3)=ds3(3)*s3(1)*s3(2)
c
c  calculate the force constants and their derivatives
c
       do i=1,3
         fdelta(i)=switch(i)*fch3
         hdelta(i)=switch(i)*hch3
         do j=1,3
           dfdelta(i,j)=dswitch(i,j)*fch3
           dhdelta(i,j)=dswitch(i,j)*hch3
         enddo
       enddo
       return
       end
c
c******************************************************
c
c
       subroutine ipforce_clnh3
c
c  calculates the symmetrised in plane bend force constants and
c  all partial derivatives involving them
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       dimension dfko(4),df1dc(4),df1dh(4)
c
c  set force constant at asymptotes
c
       f0=fkinf+ak
       f2=fkinf
c
        argfk0=bk*((rch(1)-r0ch)**2+(rch(2)-r0ch)**2
     *              +(rch(3)-r0ch)**2)
        fko=fkinf+ak*exp(-argfk0)
       fk0(1,2)=fko
       fk0(1,3)=fko
       fk0(1,4)=fko
       fk0(2,3)=fko
       fk0(2,4)=fko
       fk0(3,4)=fko
c
       do i=1,3
c
c  calc derivatives of fko wrt each of the rch(i) bonds
c
          dfko(i)=-2.0d0*ak*bk*(rch(i)-r0ch)*exp(-argfk0)
c
c  calculate the terms f1(i)
c
         arga1=aa1*rbh(i)*rbh(i)
         arga2=aa4*(rbh(i)-r0hh)*(rbh(i)-r0hh)
         a1=1.0d0-exp(-arga1)
         a2=aa2+aa3*exp(-arga2)
         f1(i)=a1*exp(-a2*(rch(i)-r0ch)**2)
c
c  and calculate the derivatives wrt rch(i) and rbh(i)
c
         duma1=2.0d0*aa1*rbh(i)*exp(-arga1)
         duma2=-2.0d0*aa3*aa4*(rbh(i)-r0hh)*exp(-arga2)
         df1dc(i)=-2.0d0*(rch(i)-r0ch)*a1*a2*exp(-a2*(rch(i)-r0ch)**2)
         df1dh(i)=duma1*exp(-a2*(rch(i)-r0ch)**2)
     *             -duma2*(rch(i)-r0ch)**2*a1*exp(-a2*(rch(i)-r0ch)**2)
       enddo
c
c  derivative of total force constant f(i,j) wrt bond length rch(k)
c  is given by dfdc(i,j,k)
c
      dfdc(1,2,1)=dfko(1)*f1(1)*f1(2)+fko*df1dc(1)*f1(2)
      dfdc(1,2,2)=dfko(2)*f1(1)*f1(2)+fko*f1(1)*df1dc(2)
      dfdc(1,2,3)=dfko(3)*f1(1)*f1(2)
c
      dfdc(1,3,1)=dfko(1)*f1(1)*f1(3)+fko*df1dc(1)*f1(3)
      dfdc(1,3,2)=dfko(2)*f1(1)*f1(3)
      dfdc(1,3,3)=dfko(3)*f1(1)*f1(3)+fko*f1(1)*df1dc(3)
c
      dfdc(2,3,1)=dfko(1)*f1(2)*f1(3)
      dfdc(2,3,2)=dfko(2)*f1(2)*f1(3)+fko*df1dc(2)*f1(3)
      dfdc(2,3,3)=dfko(3)*f1(2)*f1(3)+fko*f1(2)*df1dc(3)
c
c  derivative of total force constant f(i,j) wrt bond length rbh(k)
c  is given by dfdh(i,j,k)
c
c  only non-zero derivatives are those from rbh(i) and rbh(j)
c
       dfdh(1,2,1)=fko*df1dh(1)*f1(2)
       dfdh(1,2,2)=fko*f1(1)*df1dh(2)
       dfdh(1,2,3)=0.0d0
c
       dfdh(1,3,1)=fko*df1dh(1)*f1(3)
       dfdh(1,3,2)=0.0d0
       dfdh(1,3,3)=fko*f1(1)*df1dh(3)
c
       dfdh(2,3,1)=0.0d0
       dfdh(2,3,2)=fko*df1dh(2)*f1(3)
       dfdh(2,3,3)=fko*f1(2)*df1dh(3)
c
       return
       end
c
c******************************************************
c
c
       subroutine switchf_clnh3
c
c  calculates switching functions: s3,sphi,stheta
c  and their derivatives ds3,dsphi,dstheta
c
       implicit double precision (a-h,o-z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT1CM_clnh3/ R(N3ATOM), ENGYGS, DEGSDR(N3ATOM)
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
      COMMON/PT4CM_clnh3/ ENGYES(ISURF), DEESDR(N3ATOM,ISURF)
      COMMON/PT5CM_clnh3/ ENGYIJ(JSURF), DEIJDR(N3ATOM,JSURF)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
      COMMON/USROCM_clnh3/ PENGYGS,PENGYES(ISURF),
     +               PENGYIJ(JSURF),
     +               DGSCART(NATOM,3),DESCART(NATOM,3,ISURF),
     +               DIJCART(NATOM,3,JSURF)
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
       a1s=1.5313681d-7
       b1s=-4.6696246d0
       a2s=1.0147402d-7
       b2s=-12.363798d0
c
c  use double precision criterion:
c
c  tanh(19.0d0)=1.0d0
c
       argmax=19.0d0
       do i=1,3
         args1=a1s*(rch(i)-r0ch)*(rch(i)-b1s)**8
         if(args1.lt.argmax)then
           s1(i)=1.0d0-tanh(args1)
           ds1(i)=a1s*((rch(i)-b1s)**8
     *                 +8.0d0*(rch(i)-r0ch)*(rch(i)-b1s)**7)
           ds1(i)=-ds1(i)/cosh(args1)**2
         else
           s1(i)=0.0d0
           ds1(i)=0.0d0
         endif
c
         args2=a2s*(rch(i)-r0ch)*(rch(i)-b2s)**6
         if(args2.lt.argmax)then
           s2(i)=1.0d0-tanh(args2)
           ds2(i)=a2s*((rch(i)-b2s)**6
     *                 +6.0d0*(rch(i)-r0ch)*(rch(i)-b2s)**5)
           ds2(i)=-ds2(i)/cosh(args2)**2
         else
           s2(i)=0.0d0
           ds2(i)=0.0d0
         endif
c
c  calculate s3 and ds3
c
         args3=a3s*(rch(i)-r0ch)*(rch(i)-b3s)**2
         if (args3.lt.argmax)then
           s3(i)=1.0d0-tanh(args3)
           ds3(i)=a3s*(3.0d0*rch(i)**2-2.0d0*rch(i)*(r0ch+2.0d0*b3s)
     *          +b3s*(b3s+2.0d0*r0ch))
           ds3(i)=-ds3(i)/cosh(args3)**2
         else
           s3(i)=0.0d0
           ds3(i)=0.0d0
         endif
c
c  calculate sphi and dsphi
c
c  condition here is on the bondlength rch(i)
c  st argsphi is lt approx 19.0d0
c
         if(rch(i).lt.3.8d0)then
           argsphi=aphi*(rch(i)-r0ch)*exp(bphi*(rch(i)-cphi)**3)
           sphi(i)=1.0d0-tanh(argsphi)
           dsphi(i)=aphi*(1.0d0+3.0d0*bphi*(rch(i)-r0ch)
     *                      *(rch(i)-cphi)**2)
           dsphi(i)=dsphi(i)*exp(bphi*(rch(i)-cphi)**3)
           dsphi(i)=-dsphi(i)/cosh(argsphi)**2
         else
           sphi(i)=0.0d0
           dsphi(i)=0.0d0
         endif
c
c  calculate stheta and dstheta
c
         if(rch(i).lt.3.8d0)then
           argstheta=atheta*(rch(i)-r0ch)*exp(btheta*(rch(i)-ctheta)**3)
           stheta(i)=1.0d0-tanh(argstheta)
           dstheta(i)=atheta*(1.0d0+3.0d0*btheta*(rch(i)-r0ch)
     *           *(rch(i)-ctheta)**2)
           dstheta(i)=dstheta(i)*exp(btheta*(rch(i)-ctheta)**3)
           dstheta(i)=-dstheta(i)/cosh(argstheta)**2
         else
           stheta(i)=0.0d0
           dstheta(i)=0.0d0
         endif
       enddo
       return
       end
C
      SUBROUTINE initialize_clnh3
C
C   N3TMMN = 3 * NATOMS
C   NATOMS = the number of atoms represented by this potential function
C
C   The variable N3TMMN is the minimum value of N3TM allowed to be
C   passed by the calling routine for the number of cartesian
C   coordinates needed to represent the full system represented by this
C   potential energy surface routine.
C   N3TM must be greater than or equal to N3TMMN.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
      PARAMETER (N3TMMN = 15)
C
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
       common /ndx/     nc(3),nhb(3),nh(4,3)
C
C
C  CHECK THE NUMBER OF CARTESIAN COORDINATES SET BY THE CALLING PROGRAM
C
      IF(NATOMS.GT.25) THEN
         WRITE(NFLAG(18),1111)
 1111    FORMAT(2X,'STOP. NUMBER OF ATOMS EXCEEDS ARRAY DIMENSIONS')
         STOP
      END IF
C
       DO I=1,5
          REF(I) = ' '
       END DO
C
       REF(1)='J. Espinosa-Garcia'
       REF(2)='unpublished results, 2008'
C
      INDEXES(1) = 1
      INDEXES(2) = 7
      INDEXES(3) = 1
      INDEXES(4) = 1
      INDEXES(5) = 1
C
C
C
      IRCTNT=5
C
      CALL POTINFO_clnh3
C
      CALL ANCVRT_clnh3
c
c  calculate indexes for coordinates
c
       do ind=1,3
         icount=ind-3
         nc(ind)=3*nnc+icount
         nhb(ind)=3*nnb+icount
         do i=1,3
           nh(i,ind)=3*nnh(i)+icount
         enddo
       enddo
c
c  convert to appropriate units:
c
c  energy   in 1.0d+05 j/mol
c  time     in 1.0d-14 s
c  distance in 1.0d-10 m
c  angles   in radians
c  mass     in amu
c
       fact1=0.041840d0
       fact2=6.022045d0
c
       d1ch=d1ch*fact1
       d3ch=d3ch*fact1
       d1cb=d1cb*fact1
       d3cb=d3cb*fact1
       d1hh=d1hh*fact1
       d3hh=d3hh*fact1
       fch3=fch3*fact2
       hch3=hch3*fact2
       fkinf=fkinf*fact2
       ak=ak*fact2
C
      RETURN
      END
C
      BLOCK DATA PTPACM_clnh3
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*75 REF(5)
C
      PARAMETER(N3ATOM = 75)
      PARAMETER (ISURF = 5)
      PARAMETER (JSURF = ISURF*(ISURF+1)/2)
C
      PARAMETER (PI = 3.141592653589793D0)
      PARAMETER (NATOM = 25)
C
      COMMON/PT3CM_clnh3/ EZERO(ISURF+1)
C
      COMMON/INFOCM_clnh3/ CARTNU(NATOM,3),INDEXES(NATOM),
     +               IRCTNT,NATOMS,ICARTR,MDER,MSURF,REF
C
C
      COMMON/USRICM_clnh3/ CART(NATOM,3),ANUZERO,
     +               NULBL(NATOM),NFLAG(20),
     +               NASURF(ISURF+1,ISURF+1),NDER
C
      COMMON /POTCM_clnh3/ nnc,nnb,nnh(4),
     +               r0chr,r0chp,w1,w2,d1ch,d3ch,
     +               a1ch,b1ch,c1ch,
     +               r0hh,d1hh,d3hh,ahh,
     +               r0cb,d1cb,d3cb,acb,
     +               a3s,b3s,aphi,bphi,cphi,
     +               atheta,btheta,ctheta,
     +               fch3,hch3,
     +               fkinf,ak,bk,aa1,aa2,aa3,aa4,
     +               tau,taunh2,optau
C
       common /angles/  theta0(4,4),dtheta0(4,4,4)
       common /bonds/   rcb,rch(4),rbh(4),r0ch
       common /coords/  tcb(3),tch(4,3),tbh(4,3)
       common /delta1/  fdelta(4),hdelta(4)
       common /delta2/  dfdelta(4,4),dhdelta(4,4)
       common /force1/  fk0(4,4),f1(4),dfdc(4,4,4),dfdh(4,4,4)
       common /fsw1/    a1s,b1s,a2s,b2s
       common /ip1/     s1(4),ds1(4),s2(4),ds2(4)
       common /ndx/     nc(3),nhb(3),nh(4,3)
       common /opb1/    optheta0(4,4), doptheta0(4,4,4),ppito
       common /op1/     s3(4),ds3(4)
       common /qpdot_pl/   q(150),pdot(150)
       common /switch1/ sphi(4),dsphi(4),stheta(4),dstheta(4)
C
      DATA NASURF /1,35*0/
      DATA NDER /1/
       DATA NFLAG /1,1,15*0,6,0,0/
C
      DATA ANUZERO /0.0D0/
      DATA ICARTR,MSURF,MDER/1,0,1/
      DATA NULBL /25*0/
      DATA NATOMS /5/
C
       DATA nnc     /2/
       DATA nnb     /5/
       DATA nnh     /3,4,1,0/
       DATA r0chr   /   1.01410d0/
       DATA r0chp   /   1.02700d0/
       DATA w1      /   1.00000d0/
       DATA w2      /   1.01400d0/
       DATA d1ch    / 119.058d0/
       DATA d3ch    /  20.000d0/
       DATA a1ch    /   2.125000d0/
       DATA b1ch    /  -0.090000d0/
       DATA c1ch    /   22.00000d0/
       DATA r0hh    /   1.27730d0/
       DATA d1hh    / 109.850d0/
       DATA d3hh    /  18.400d0/
       DATA ahh     /   1.8600d0/
       DATA r0cb    /   2.10400d0/
       DATA d1cb    /  65.100d0/
       DATA d3cb    /  16.530d0/
       DATA acb     /   0.7780000d0/
       DATA a3s     /   1.0897000d0/
       DATA b3s     /  -0.8088000d0/
       DATA aphi    /   6.7730500d0/
       DATA bphi    /   6.8000000d0/
       DATA cphi    /   1.9226100d0/
       DATA atheta  /   6.7359700d0/
       DATA btheta  /   6.7000000d0/
       DATA ctheta  /   1.9505500d0/
       DATA fch3    /  -0.0520000d0/
       DATA hch3    /  -0.0510000d0/
       DATA fkinf   /   0.6950000d0/
       DATA ak      /  -0.0100000d0/
       DATA bk      /   0.1000100d0/
       DATA aa1     /   3.503370d0/
       DATA aa2     /   6.130490d0/
       DATA aa3     /   6.100000d0/
       DATA aa4     /   3.232430d0/
       DATA tau     /   1.9022600d0/
       DATA taunh2  /   1.8046700d0/
       DATA optau   /   1.1847700d0/
C
       END

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


! { dg-do run }

  integer, allocatable :: a, b(:), c(:,:)
  logical :: l
  if (allocated (a) .or. allocated (b) .or. allocated (c)) call abort

!$omp parallel private (a, b, c, l)
  l = .false.
  if (allocated (a) .or. allocated (b) .or. allocated (c)) call abort

!$omp single
  allocate (a, b(6:9), c(3, 8:9))
  a = 4
  b = 5
  c = 6
!$omp end single copyprivate (a, b, c)

  if (.not.allocated (a)) call abort
  if (.not.allocated (b) .or. size (b) /= 4) call abort
  if (lbound (b, 1) /= 6 .or. ubound (b, 1) /= 9) call abort
  if (.not.allocated (c) .or. size (c) /= 6) call abort
  if (size (c, 1) /= 3 .or. size (c, 2) /= 2) call abort
  if (lbound (c, 1) /= 1 .or. ubound (c, 1) /= 3) call abort
  if (lbound (c, 2) /= 8 .or. ubound (c, 2) /= 9) call abort
  if (a /= 4 .or. any (b /= 5) .or. any (c /= 6)) call abort

!$omp single
  deallocate (a, b, c)
  if (allocated (a) .or. allocated (b) .or. allocated (c)) call abort
  allocate (a, b(0:4), c(3, 2:7))
  a = 1
  b = 2
  c = 3
!$omp end single copyprivate (a, b, c)

  if (.not.allocated (a)) call abort
  if (.not.allocated (b) .or. size (b) /= 5) call abort
  if (lbound (b, 1) /= 0 .or. ubound (b, 1) /= 4) call abort
  if (.not.allocated (c) .or. size (c) /= 18) call abort
  if (size (c, 1) /= 3 .or. size (c, 2) /= 6) call abort
  if (lbound (c, 1) /= 1 .or. ubound (c, 1) /= 3) call abort
  if (lbound (c, 2) /= 2 .or. ubound (c, 2) /= 7) call abort
  if (a /= 1 .or. any (b /= 2) .or. any (c /= 3)) call abort

!$omp single
  l = .true.
  deallocate (a, b, c)
  if (allocated (a) .or. allocated (b) .or. allocated (c)) call abort
  allocate (a, b(2:6), c(3:5, 3:8))
  a = 7
  b = 8
  c = 9
!$omp end single copyprivate (a, b, c)

  if (.not.allocated (a)) call abort
  if (.not.allocated (b) .or. size (b) /= 5) call abort
  if (l) then
    if (lbound (b, 1) /= 2 .or. ubound (b, 1) /= 6) call abort
  else
    if (lbound (b, 1) /= 0 .or. ubound (b, 1) /= 4) call abort
  end if
  if (.not.allocated (c) .or. size (c) /= 18) call abort
  if (size (c, 1) /= 3 .or. size (c, 2) /= 6) call abort
  if (l) then
    if (lbound (c, 1) /= 3 .or. ubound (c, 1) /= 5) call abort
    if (lbound (c, 2) /= 3 .or. ubound (c, 2) /= 8) call abort
  else
    if (lbound (c, 1) /= 1 .or. ubound (c, 1) /= 3) call abort
    if (lbound (c, 2) /= 2 .or. ubound (c, 2) /= 7) call abort
  end if
  if (a /= 7 .or. any (b /= 8) .or. any (c /= 9)) call abort

!$omp end parallel
end


text:
! { dg-do run }
! PR47285 G format outputs wrong number of characters.
! Test case prepared by Jerry DeLisle <jvdelisle@gcc.gnu.org>
       PROGRAM FOO
       character(len=50) :: buffer

       WRITE(buffer,"(G0.5,'<')") -10000.
       if (buffer.ne."-10000.<") call abort
       WRITE(buffer,"(G1.5E5,'<')") -10000.
       if (buffer.ne."*<") call abort
       WRITE(buffer,"(G2.5E5,'<')") -10000.
       if (buffer.ne."**<") call abort
       WRITE(buffer,"(G3.5E5,'<')") -10000.
       if (buffer.ne."***<") call abort
       WRITE(buffer,"(G4.5E5,'<')") -10000.
       if (buffer.ne."****<") call abort
       WRITE(buffer,"(G5.5E5,'<')") -10000.
       if (buffer.ne."*****<") call abort
       WRITE(buffer,"(G6.5E5,'<')") -10000.
       if (buffer.ne."******<") call abort
       WRITE(buffer,"(G7.5E5,'<')") -10000.
       if (buffer.ne."*******<") call abort
       WRITE(buffer,"(G8.5E5,'<')") -10000.
       if (buffer.ne."********<") call abort
       WRITE(buffer,"(G9.5E5,'<')") -10000.
       if (buffer.ne."*********<") call abort
       WRITE(buffer,"(G10.5E5,'<')") -10000.
       if (buffer.ne."**********<") call abort
       WRITE(buffer,"(G11.5E5,'<')") -10000.
       if (buffer.ne."***********<") call abort
       WRITE(buffer,"(G12.5E5,'<')") -10000.
       if (buffer.ne."************<") call abort
       WRITE(buffer,"(G13.5E5,'<')") -10000.
       if (buffer.ne."*************<") call abort
       WRITE(buffer,"(G14.5E5,'<')") -10000.
       if (buffer.ne."-10000.       <") call abort
       WRITE(buffer,"(G15.5E5,'<')") -10000.
       if (buffer.ne." -10000.       <") call abort
       WRITE(buffer,"(G16.5E5,'<')") -10000.
       if (buffer.ne."  -10000.       <") call abort

       STOP
       END

ext: f


text:
! signal.f90
!
! Example program that registers a signal handler for SIGINT.
!
! Author:  Philipp Engel
! Licence: ISC
program main
    use, intrinsic :: iso_c_binding
    use :: unix
    implicit none
    integer        :: rc
    type(c_funptr) :: ptr

    ! Register signal handler.
    ptr = c_signal(SIGINT, c_funloc(sigint_handler))

    print '(a)', 'Press CTRL + C to send SIGINT.'

    do
        print '(a)', 'zzz ...'
        rc = c_usleep(10**6)
    end do
contains
    subroutine sigint_handler(signum) bind(c)
        !! Signal handler for SIGINT.
        integer(kind=c_int), intent(in), value :: signum

        print '(a, i0, a)', 'Received SIGINT (', signum, '). Terminating ...'
        stop
    end subroutine sigint_handler
end program main

max_stars_repo_name: interkosmos/fortran-unix

max_issues_repo_name: interkosmos/fortran-unix

max_forks_repo_name: interkosmos/fortran-unix


text:
MODULE m_tetsrt

   IMPLICIT NONE

   CONTAINS

   PURE FUNCTION tetsrt(etetra) Result(ind)

      REAL,             INTENT(IN)     :: etetra(:)

      INTEGER :: ind(SIZE(etetra))

      INTEGER i,j
      INTEGER tmp

      DO i = 1, SIZE(etetra)
         ind(i) = i
      ENDDO

      !Sort the energies in the tetrahedron in ascending order
      DO i = 1, SIZE(etetra)-1
         DO j = i+1, SIZE(etetra)
            IF (etetra(ind(i)).GT.etetra(ind(j))) THEN
               tmp = ind(i)
               ind(i) = ind(j)
               ind(j) = tmp
            ENDIF
         ENDDO
      ENDDO

   END FUNCTION tetsrt

END MODULE m_tetsrt

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
 A photographer

JulioD is a guy that lived in a very big, noisy and contaminated city for over 30 years. As destiny will have it, he then moved to Davis. This a photographic chronicle of his live in Davis, Ca. I will also try to start a night photography section.

Day job

Currently, I am doing my last postdoc at UC Davis. My area of research is structural biology with a lot of electron microscopy and 2d crystallography. This means that I spent a good part of my day taking pictures (or rather projections) of things with sizes of about 1 billionth of a meter. But it is fun. Althought, the microscopes are noisy but with a good set of earphones that is no problem at all.

Previously, I worked doing wet biochemistry and some theoretical stuff, mostly molecular dynamics simulations and docking predictions. I actually created a few videos (examples: http://www.youtube.com/watch?vo5a39tbT9w or http://www.youtube.com/watch?vQ0VEbNaCPM).

Speaking about noise, it is rather unusual for me to be seen without earphones (or headphones). I usually wear set Bose Q3s http://www.headphoneinfo.com/content/BostQuietComfort3HeadphonesReview.htm which are awesome for noise environments (chatting people, rotary pumps, A/C, airplane travel). Downside of such headphones is they run on a proprietary battery. Thats why I got a pair Shure SE530 http://www.headphoneinfo.com/content/ShureSE530HeadphonesReview825.htm. No, I am not trying to sell earphones (or iPods, for that matter). And no, I am no audiophile. It is just that after a certain level (~100 bucks) good (not audiophile good but good enough for my purposes) options are limited. But enough about that.

I used to have a Cybershot (up until december 2009). Then, I moved to a beautiful DSLR. So, if you see any shots taken before december 2009, the Cybershot did it. Afterwards, it is all Nikons fault.

What I like about Davis is...

The night sky is beautiful. For example, last night (06/27/2010) the moon looked like this:

http://lh4.ggpht.com/_fP4_LG3LE9Q/TCe9s5kpvoI/AAAAAAAACSw/nnVcLXkeEDg/s800/DSC_4219.jpg

Info: ISO 640, 170 mm, f/5.6 at 1/1250 sec.

 The photographer production

Fireworks and a patriotic cow? (07/04/2010)

>Image(fireworks.jpg, thumbnail, 300, ISO 200, 35 mm, f/22 and 3.2 seconds. Manual exposure and focus)<

This is the first picture I shot in Davis, Ca on 12/29/2007:

>Image(number1.jpg, thumbnail, 300, This was taken from a window of the now closed Davis Inn.)<

Night photography

Last night (06282010) Scorpio was above the horizont and I capture a fairly recognizable part of it:

This month Ursa Major has been high in the sky. http://en.wikipedia.org/wiki/Ursa_Major This is a failed attempt to capture it (06/16/2010):

>Image(Attempt1.jpg, left, thumbnail, 300, ISO 200, 55 mm, f/5.6 and 10.0 seconds)<

If you care about photographic info this pic was taken with XISO 200, 200 mm, f/5.6 at 1/30 secX  (corrected, thanks for pointing out the mistake wl).

The night of june 28th 2010 I shot this:

>Image(UrsaMajor.jpg, left, thumbnail, 400, Here we go, the brightest part of Ursa Major. ISO 200, 18 mm, f/4.5 and 15 seconds. Manual focus.)<

You dont have any shots of the pleiades do you? Daubert

Looking through my library it seems that back in april 16th 2010 I shot this:

>Image(Moonand7sis.jpg, thumbnail, 500, The moon and the seven sisters (the pleiades). ISO 200, 200 mm, f/5.6 and 2.5 seconds. Auto focus, I think.)<

20100629 13:50:56 nbsp I love that last photo of the crescent moon.  I am anxiously awaiting more night photo shots.  A great place for night photography in Davis is the Davis cemetery.  Very little light pollution, few people, and a very peaceful place to be. Users/DagonJones

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! NHDS is a numerical solver to the linear hot-plasma dispersion relation
! Copyright (C) 2018 Daniel Verscharen (d.verscharen@ucl.ac.uk)
!All rights reserved.
!
!Redistribution and use in source and binary forms, with or without
!modification, are permitted provided that the following conditions are met:
!
!1. Redistributions of source code must retain the above copyright notice, this
!   list of conditions and the following disclaimer.
!2. Redistributions in binary form must reproduce the above copyright notice,
!   this list of conditions and the following disclaimer in the documentation
!   and/or other materials provided with the distribution.
!
!THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
!ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
!WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
!DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
!ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
!(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
!LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
!ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
!(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!The views and conclusions contained in the software and documentation are those
!of the authors and should not be interpreted as representing official policies,
!either expressed or implied, of the NHDS project.

program NHDS
use input_params
implicit none
integer :: j,i,l
double complex :: x
double precision, allocatable, dimension(:) :: kk,theta
double complex, allocatable, dimension(:) :: guesses
double precision :: dk, dth
character*30 :: filename

dk=0.; dth=0.
call get_command_argument(1,filename); filename=adjustl(filename)
call set_parameters(filename)

allocate(kk(ksteps))
if (kth_file) then
   ! if reading k,th values from file
   ! ksteps=theta_steps=number of (k,th) pairs in input file
   allocate(theta(ksteps),guesses(ksteps))
   kth_filename=adjustl(kth_filename)
   open(unit=793,file=kth_filename)
   read(793,*) kk
   read(793,*) theta
   read(793,*) guesses
   close(793)
else
   allocate(theta(theta_steps))
   ! If using range of k and theta, create the kk, theta arrays to use
   if (ksteps .gt. 1) dk = (krange(2)-krange(1))/(1.d0*ksteps-1.d0)
   if (theta_steps .gt. 1) dth = (theta_range(2)-theta_range(1))/(1.d0*theta_steps-1.d0)
   do i=1,ksteps
      kk(i)=krange(1)+dk*(1.d0*i-1.d0)
    enddo
   do i=1,theta_steps
      theta(i)=theta_range(1)+dth*(1.d0*i-1.d0)
   enddo
endif

do j=1,numspec
  Omega(j)=charge(j)/mass(j)
  ell(j)=sqrt(mass(j)/(density(j)*charge(j)*charge(j)))
  vtherm(j)=sqrt(beta(j)/(density(j)*mass(j)))
enddo

open(unit=10,file='output_'//trim(filename)//'.dat',status='replace',action='write')
write (10,*) '# kk, theta, omega, gamma, Re(Ey/Ex), Im(Ey/Ex), Re(Ez/Ex), Im(Ez/Ex), energy, quality'

if (output_mom) then
   open(unit=11,file='Plasma_'//trim(filename)//'.dat',status='replace',action='write')
   write (11,*) '# kk, theta, Re(dn), Im(dn), Re(dUperpx), Im(dUperpx), Re(dUperpy), Im(dUperpy), Re(dUpar), Im(dUpar)'
endif
if (output_EB) then
   open(unit=12,file='EB_'//trim(filename)//'.dat',status='replace',action='write')
   write (12,*) '# kk, theta, Re(Ex), Im(Ex), Re(Ey), Im(Ey), Re(Ez), Im(Ez), Re(Bx), Im(Bx), Re(By), Im(By), Re(Bz), Im(Bz)'
endif

if (kth_file) then
   do i=1,ksteps
      x = guesses(i)
      call compute(kk(i),theta(i),x,output_mom,output_EB)
   enddo
else
   do l=1,theta_steps
      x=initial_guess ! guess w in units of wci:
      do i=1,ksteps
         call compute(kk(i),theta(l),x, output_mom,output_EB)
         if (i .eq. 1) initial_guess=x
      enddo
   enddo
endif

close(10)
close(11)
close(12)

end program

subroutine compute(kk,theta,x,outputm,outputeb)
   use input_params
   double precision, intent(in) :: kk, theta
   double complex, intent(inout) :: x
   logical, intent(in) :: outputm, outputeb
   double complex :: pol,polz,xi,deltaRj, Ek(3), Bk(3)
   double complex :: dUperpx,dUperpy,dUpar, Avec(3)
   double precision :: kperp,kz,energy,gamma_contribution(10)
   real :: quality
   integer :: j
   kz=kk*cos(theta*M_PI/180.d0)
   kperp=kk*sin(theta*M_PI/180.d0)

   call newton_method(x,kz,kperp,quality)
   call calc_polarization(pol,polz,x,kz,kperp)
   call waveenergy(energy,x,kz,kperp,gamma_contribution)

!  if (i .eq. 1) initial_guess=x

   if (outputm) then
     do j=1,numspec
       call calc_xi(xi,j,pol,polz,x,kz,kperp)
       ! second parameter is index of species
       call calc_fluctRj(deltaRj,dUperpx,dUperpy,dUpar,j,pol,polz,x,kz,kperp)
       ! fifth parameter is index of species
       write(11,'(10F14.5)') kk, theta, real(xi), aimag(xi), real(dUperpx), &
                            aimag(dUperpx), real(dUperpy), aimag(dUperpy), &
                            real(dUpar), aimag(dUpar)
     enddo
     write(11,*) ''
   endif

   if (outputeb) then

     if (ampl_mode.EQ.1) then
     	! Avec fulfills the equation  c * dE = Avec * vA * dB_x
     	Avec(1)=uniti*x/(kz*pol)
     	Avec(2)=-x/kz
     	Avec(3)=polz*uniti*x/(kz*pol)
     else if (ampl_mode.EQ.2) then
     	! Avec fulfills the equation  c * dE = Avec * vA * dB_y
     	Avec(1)=x/(kz-polz*kperp)
     	Avec(2)=uniti*pol*x/(kz-polz*kperp)
     	Avec(3)=x/((kz/polz)-kperp)
     else if (ampl_mode.EQ.3) then
     	! Avec fulfills the equation  c * dE = Avec * vA * dB_z
     	Avec(1)=-uniti*x/(kperp*pol)
     	Avec(2)=x/kperp
     	Avec(3)=-uniti*x*polz/(kperp*pol)
     endif

      Ek(1)=Avec(1)*vAc*ampl
      Ek(2)=Avec(2)*vAc*ampl
      Ek(3)=Avec(3)*vAc*ampl

      Bk(1)=-kz*Ek(2)/(x*vAc)
      Bk(2)=(kz*Ek(1)-kperp*Ek(3))/(x*vAc)
      Bk(3)=kperp*Ek(2)/(x*vAc)

      write(12,'(14E14.5)') kk, theta,   real(Ek(1)), aimag(Ek(1)), real(Ek(2))&
                          ,aimag(Ek(2)), real(Ek(3)), aimag(Ek(3)), real(Bk(1))&
                          ,aimag(Bk(1)), real(Bk(2)), aimag(Bk(2)), real(Bk(3))&
                          ,aimag(Bk(3))
   endif

!  write (*, '(10F9.5)') kk(i),theta(j),real(x),aimag(x),real(pol),aimag(pol),&
!                        real(polz),aimag(polz),energy,quality
   write (*, '(3F10.5,3E14.5)' ) kk, theta, real(x), aimag(x), quality
   write (10,'(3F10.5,7E14.5)') kk, theta, real(x), aimag(x), real(pol*uniti), aimag(pol*uniti),&
                          real(polz), aimag(polz), energy, quality

! if (i.EQ.10) call write_delta_f(1,kz,x,pol,polz) ! first parameter is index of species
end subroutine compute

max_stars_repo_name: tulasinandan/NHDS

max_issues_repo_name: tulasinandan/NHDS

max_forks_repo_name: tulasinandan/NHDS


text:
      program fkinDiagon_kry_p
c     ----------------------------------------------------------------
c     $Revision: 1.2 $
c     $Date: 2009/09/30 23:42:12 $
c     ----------------------------------------------------------------
c     Programmer(s): Allan G. Taylor, Alan C. Hindmarsh and
c                    Radu Serban @ LLNL
c     ----------------------------------------------------------------
c     Simple diagonal test with Fortran interface, using
c     user-supplied preconditioner setup and solve routines (supplied
c     in Fortran, below).
c
c     This example does a basic test of the solver by solving the
c     system:
c               f(u) = 0  for
c               f(u) = u(i)^2 - i^2
c
c      No scaling is done.
c      An approximate diagonal preconditioner is used.
c
c      Execution command: mpirun -np 4 fkinDiagon_kry_p
c     ----------------------------------------------------------------
c
      implicit none

      include "mpif.h"

      integer ier, size, globalstrat, rank, mype, npes
      integer maxl, maxlrst
      integer*4 localsize
      parameter(localsize=32)
      integer*4 neq, nlocal, msbpre, baseadd, i, ii
      integer*4 iout(15)
      double precision rout(2)
      double precision pp, fnormtol, scsteptol
      double precision uu(localsize), scale(localsize)
      double precision constr(localsize)

      common /pcom/ pp(localsize), mype, npes, baseadd, nlocal

      nlocal = localsize
      neq = 4 * nlocal
      globalstrat = 0
      fnormtol = 1.0d-5
      scsteptol = 1.0d-4
      maxl = 10
      maxlrst = 2
      msbpre  = 5

c     The user MUST call mpi_init, Fortran binding, for the fkinsol package
c     to work. The communicator, MPI_COMM_WORLD, is the only one common
c     between the Fortran and C bindings. So in the following, the communicator
c     MPI_COMM_WORLD is used in calls to mpi_comm_size and mpi_comm_rank
c     to determine the total number of processors and the rank (0 ... size-1)
c     number of this process.

      call mpi_init(ier)
      if (ier .ne. 0) then
         write(6,1210) ier
 1210    format('MPI_ERROR: MPI_INIT returned IER = ', i2)
         stop
      endif

      call fnvinitp(mpi_comm_world, 3, nlocal, neq, ier)
      if (ier .ne. 0) then
         write(6,1220) ier
 1220    format('SUNDIALS_ERROR: FNVINITP returned IER = ', i2)
         call mpi_finalize(ier)
         stop
      endif

      call mpi_comm_size(mpi_comm_world, size, ier)
      if (ier .ne. 0) then
         write(6,1222) ier
 1222    format('MPI_ERROR: MPI_COMM_SIZE returned IER = ', i2)
         call mpi_abort(mpi_comm_world, 1, ier)
         stop
      endif

      if (size .ne. 4) then
         write(6,1230)
 1230    format('MPI_ERROR: must use 4 processes')
         call mpi_finalize(ier)
         stop
      endif
      npes = size

      call mpi_comm_rank(mpi_comm_world, rank, ier)
      if (ier .ne. 0) then
         write(6,1224) ier
 1224    format('MPI_ERROR: MPI_COMM_RANK returned IER = ', i2)
         call mpi_abort(mpi_comm_world, 1, ier)
         stop
      endif

      mype = rank
      baseadd = mype * nlocal

      do 20 ii = 1, nlocal
         i = ii + baseadd
         uu(ii) = 2.0d0 * i
         scale(ii) = 1.0d0
         constr(ii) = 0.0d0
 20   continue

      call fkinmalloc(iout, rout, ier)

      if (ier .ne. 0) then
         write(6,1231)ier
 1231    format('SUNDIALS_ERROR: FKINMALLOC returned IER = ', i2)
         call mpi_abort(mpi_comm_world, 1, ier)
         stop
      endif

      call fkinsetiin('MAX_SETUPS', msbpre, ier)
      call fkinsetrin('FNORM_TOL', fnormtol, ier)
      call fkinsetrin('SSTEP_TOL', scsteptol, ier)
      call fkinsetvin('CONSTR_VEC', constr, ier)

      call fkinspgmr(maxl, maxlrst, ier)
      call fkinspilssetprec(1, ier)

      if (mype .eq. 0) write(6,1240)
 1240 format('Example program fkinDiagon_kry_p:'//
     1       ' This FKINSOL example',
     2       ' solves a 128 eqn diagonal algebraic system.'/
     3       ' Its purpose is to demonstrate the use of the Fortran',
     4       ' interface'/' in a parallel environment.')

      call fkinsol(uu, globalstrat, scale, scale, ier)
      if (ier .lt. 0) then
         write(6,1242) ier, iout(9)
 1242    format('SUNDIALS_ERROR: FKINSOL returned IER = ', i2, /,
     1          '                Linear Solver returned IER = ', i2)
         call mpi_abort(mpi_comm_world, 1, ier)
         stop
      endif

      if (mype .eq. 0) write(6,1245) ier
 1245 format(/' FKINSOL return code is ', i4)

      if (mype .eq. 0) write(6,1246)
 1246 format(/' The resultant values of uu (process 0) are:'/)

      do 30 i = 1, nlocal, 4
         if(mype .eq. 0) write(6,1256) i + baseadd, uu(i), uu(i+1),
     1                                 uu(i+2), uu(i+3)
 1256    format(i4, 4(1x, f10.6))
 30   continue

      if (mype .eq. 0) write(6,1267) iout(3), iout(14), iout(4),
     1                               iout(12), iout(13), iout(15)
 1267 format(/'Final statistics:'//
     1     '  nni = ', i3, ',  nli  = ', i3, /,
     2     '  nfe = ', i3, ',  npe  = ', i3, /,
     3     '  nps = ', i3, ',  ncfl = ', i3)

      call fkinfree

c     An explicit call to mpi_finalize (Fortran binding) is required by
c     the constructs used in fkinsol.
      call mpi_finalize(ier)

      stop
      end

c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     The function defining the system f(u) = 0 must be defined by a Fortran
c     function with the following name and form.

      subroutine fkfun(uu, fval, ier)

      implicit none

      integer mype, npes, ier
      integer*4 baseadd, nlocal, i, localsize
      parameter(localsize=32)
      double precision pp
      double precision fval(*), uu(*)

      common /pcom/ pp(localsize), mype, npes, baseadd, nlocal

      do 10 i = 1, nlocal
 10      fval(i) = uu(i) * uu(i) - (i + baseadd) * (i + baseadd)

      return
      end

c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     The routine kpreco is the preconditioner setup routine. It must have
c     that specific name be used in order that the c code can find and link
c     to it.  The argument list must also be as illustrated below:

      subroutine fkpset(udata, uscale, fdata, fscale,
     1                  vtemp1, vtemp2, ier)

      implicit none

      integer ier, mype, npes
      integer*4 localsize
      parameter(localsize=32)
      integer*4 baseadd, nlocal, i
      double precision pp
      double precision udata(*), uscale(*), fdata(*), fscale(*)
      double precision vtemp1(*), vtemp2(*)

      common /pcom/ pp(localsize), mype, npes, baseadd, nlocal

      do 10 i = 1, nlocal
 10      pp(i) = 0.5d0 / (udata(i)+ 5.0d0)

      ier = 0

      return
      end

c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     The routine kpsol is the preconditioner solve routine. It must have
c     that specific name be used in order that the c code can find and link
c     to it.  The argument list must also be as illustrated below:

      subroutine fkpsol(udata, uscale, fdata, fscale,
     1                  vv, ftem, ier)

      implicit none

      integer ier, mype, npes
      integer*4 baseadd, nlocal, i
      integer*4 localsize
      parameter(localsize=32)
      double precision udata(*), uscale(*), fdata(*), fscale(*)
      double precision vv(*), ftem(*)
      double precision pp

      common /pcom/ pp(localsize), mype, npes, baseadd, nlocal

      do 10 i = 1, nlocal
 10      vv(i) = vv(i) * pp(i)

      ier = 0

      return
      end

max_stars_repo_name: gregmedlock/roadrunnerwork

max_issues_repo_name: gregmedlock/roadrunnerwork

max_forks_repo_name: gregmedlock/roadrunnerwork


text:
c
c     version 1.2  Mar 1992
c
c     ----------------------------------------------------------------
c     |                                                              |
c     |              ********   R E D U C E  ********                |
c     |                                                              |
c     |            Change the reducing velocity applied              |
c     |                      to a "tx.in" file                       |
c     |                                                              |
c     |                   Written by C. A. Zelt                      |
c     |                                                              |
c     |                Geological Survey of Canada                   |
c     |                  Ottawa, Canada K1A 0Y3                      |
c     |                                                              |
c     ----------------------------------------------------------------
c
c     Note: the old or new reducing velocity may be zero for unreduced
c           data
c
      open(unit=10, file='tx.in')
      open(unit=11, file='tx.out')
c
      write(6,55)
55    format(/'Enter old and new reducing velocities (km/s)')
      read(5,*) vred1,vred2
c
      if(vred1.eq.0.) then
        rvred1=0.
      else
        rvred1=1./vred1
      end if
      if(vred2.eq.0.) then
        rvred2=0.
      else
        rvred2=1./vred2
      end if
c
100   read(10,*,end=999) xf,tf,uf,if
5     format(3f10.3,i10)
c
      if(if.eq.0) xshot=xf
      if(if.gt.0) tf=tf+abs(xf-xshot)*(rvred1-rvred2)
      write(11,5) xf,tf,uf,if
      if(if.eq.-1) go to 999
c
      go to 100
c
999   stop
      end