text:
      module dihedral_module

c***********************************************************************
c
c     dl_poly module for defining dihedral potential arrays
c     copyright - daresbury laboratory
c
c     author  - w. smith     sep 2003
c     adapted - p.-a. cazade oct 2007 : solvation, free energy, etc.
c     adapted - w.smith      jan 2011 : metadynamics
c
c***********************************************************************

      use config_module
      use error_module
      use metafreeze_module
      use parse_module
      use property_module
      use setup_module
      use site_module
      use solvation_module
      use utility_module
      use vdw_module

      implicit none

      real(8), allocatable :: prmdih(:,:)
      integer, allocatable :: listdih(:,:)
      integer, allocatable :: numdih(:),keydih(:),lstdih(:,:)

      save prmdih,listdih,numdih,keydih,lstdih

      contains

      subroutine alloc_dih_arrays(idnode,mxnode)

      implicit none

      logical safe
      integer i,fail,idnode,mxnode
      dimension fail(5)

      safe=.true.

c     allocate arrays

      fail(:)=0

      allocate (prmdih(mxtdih,mxpdih),stat=fail(1))
      allocate (numdih(mxtmls),stat=fail(2))
      allocate (keydih(mxtdih),stat=fail(3))
      allocate (lstdih(mxtdih,4),stat=fail(4))
      allocate (listdih(mxdihd,5),stat=fail(5))

      if(any(fail.gt.0))safe=.false.
      if(mxnode.gt.1)call gstate(safe)
      if(.not.safe)call error(idnode,1011)

c     initialise numdih array

      do i=1,mxtmls
        numdih(i)=0
      enddo

      end subroutine alloc_dih_arrays

      subroutine define_dihedrals
     x  (safe,idnode,itmols,ndihed,nsite,engunit)

c***********************************************************************
c
c     dl_poly subroutine for defining dihedral angles
c
c     copyright - daresbury laboratory
c     author    - w. smith    aug 2003
c
c***********************************************************************

      implicit none

      logical safe
      character*8 keyword
      character*1 message(80)
      integer idnode,itmols,ndihed,nsite,ntmp,idih,j,keytmp
      integer iatm1,iatm2,iatm3,iatm4,idum,isite1,isite2,isite3
      integer isite4
      real(8) engunit,parpot(mxpdih)

      ntmp=intstr(record,lenrec,idum)
      numdih(itmols)=numdih(itmols)+ntmp

      if(idnode.eq.0)then
        write(nrite,"(/,1x,'number of dihedral angles',
     x    6x,i10)")numdih(itmols)
        write(nrite,"(/,1x,'dihedral angle details:',
     x    /,/,12x,'unit',5x,'key',5x,'index',5x,'index',5x,'index',
     x    5x,'index',7x,'f-const',8x,'angle',9x,'trig',11x,
     x    '1-4 elec',7x,'1-4 vdw',/)")
      endif

      do idih=1,ntmp

        ndihed=ndihed+1
        if(ndihed.gt.mxtdih)call error(idnode,60)

c     read dihedral bond angle potential parameters

        call getrec(safe,idnode,nfield)
        if(.not.safe)return

        call copystring(record,message,80)
        call lowcase(record,4)
        call getword(keyword,record,4,lenrec)

        if(keyword(1:4).eq.'cos ')then
          keytmp=1
        elseif(keyword(1:4).eq.'harm')then
          keytmp=2
        elseif(keyword(1:4).eq.'hcos')then
          keytmp=3
        elseif(keyword(1:4).eq.'cos3')then
          keytmp=4
        elseif(keyword(1:4).eq.'ryck')then
          keytmp=5
        elseif(keyword(1:4).eq.'rbf')then
          keytmp=6
        elseif(keyword(1:4).eq.'opls')then
          keytmp=7
        else
          if(idnode.eq.0)write(nrite,*)message
          call error(idnode,448)
        endif

        iatm1=intstr(record,lenrec,idum)
        iatm2=intstr(record,lenrec,idum)
        iatm3=intstr(record,lenrec,idum)
        iatm4=intstr(record,lenrec,idum)
        parpot(1)=dblstr(record,lenrec,idum)
        parpot(2)=dblstr(record,lenrec,idum)
        parpot(3)=dblstr(record,lenrec,idum)
        parpot(4)=dblstr(record,lenrec,idum)
        parpot(5)=dblstr(record,lenrec,idum)

        isite1=nsite-numsit(itmols)+iatm1
        isite2=nsite-numsit(itmols)+iatm2
        isite3=nsite-numsit(itmols)+iatm3
        isite4=nsite-numsit(itmols)+iatm4

c     test for frozen atom pairs

        if(lfzsit(isite1)*lfzsit(isite2)*
     x    lfzsit(isite3)*lfzsit(isite4).ne.0)then

          numdih(itmols)=numdih(itmols)-1
          if(idnode.eq.0)
     x      write(nrite,"(4x,a8,i10,4x,a4,4i10,10f15.6)")
     x      '*frozen*',idih,keyword(1:4),iatm1,iatm2,iatm3,iatm4,
     x      (parpot(j),j=1,mxpdih)

        else

          if(idnode.eq.0)
     x      write(nrite,"(12x,i10,4x,a4,4i10,10f15.6)")
     x      idih,keyword(1:4),iatm1,iatm2,iatm3,iatm4,
     x      (parpot(j),j=1,mxpdih)

        endif

        keydih(ndihed)=keytmp
        lstdih(ndihed,1)=iatm1
        lstdih(ndihed,2)=iatm2
        lstdih(ndihed,3)=iatm3
        lstdih(ndihed,4)=iatm4
        prmdih(ndihed,1)=parpot(1)
        prmdih(ndihed,2)=parpot(2)
        prmdih(ndihed,3)=parpot(3)
        prmdih(ndihed,4)=parpot(4)
        prmdih(ndihed,5)=parpot(5)

c     convert energies to internal units and angles to radians

        prmdih(ndihed,1)=prmdih(ndihed,1)*engunit

        if(keytmp.eq.4)then

          prmdih(ndihed,2)=prmdih(ndihed,2)*engunit
          prmdih(ndihed,3)=prmdih(ndihed,3)*engunit

        elseif(keytmp.eq.7)then

          prmdih(ndihed,2)=prmdih(ndihed,2)*engunit
          prmdih(ndihed,3)=prmdih(ndihed,3)*engunit
          prmdih(ndihed,4)=prmdih(ndihed,4)*engunit
          prmdih(ndihed,5)=prmdih(ndihed,5)*(pi/180.d0)

        else

          prmdih(ndihed,2)=prmdih(ndihed,2)*(pi/180.d0)

        endif

      enddo

      return
      end subroutine define_dihedrals

      subroutine dihfrc
     x  (lsolva,lfree,lexcite,idnode,imcon,mxnode,ntdihd,keyfce,
     x  dlrpot,epsq,engcpe,engdih,engsrp,rcut,rvdw,alpha,vircpe,
     x  virdih,virsrp)

c***********************************************************************
c
c     dl_poly subroutine for calculating dihedral energy and force
c     terms in molecular dynamics.
c
c     version 3: scale factors for reduces electrostatic and vdw
c     1-4 interactions.
c
c     NOTE: assumes 1-4 interactions are in the exclude list
c
c     block as opposed to stride version
c
c     copyright - daresbury laboratory
c     author    - w. smith       mar 1992
c     modified  - t. forester    dec 1993
c     modified  - t. forester    jun 1995 - stress tensor added
c     modified  - a. smondyrev   may 2000 - ryckaert-bellemans potentials
c     modified  - p.-a. cazade   oct 2007 - solvation etc
c
c***********************************************************************

      implicit none

      logical safe,lsolva,lfree,lexcite,lselect
      logical idrive,jdrive,kdrive,ldrive
      integer i,k,ii,kk,ntdihd,idnode,mxnode,idih1,idih2,ia,ib
      integer ic,id,imcon,ka,kb,l,keyfce,fail1,fail2,fail3,kkk
      real(8) phi,twopi,rtwopi,dterm,srpot
      real(8) engdih,virdih,engc14,engs14,virs14,rrbc,xab,yab,erc,fer
      real(8) zab,xbc,ybc,zbc,xcd,ycd,zcd,pbx,pby,pbz,pb2,rpb1,rpb2
      real(8) pcx,pcy,pcz,pc2,rpc1,rpc2,pbpc,cosp,sinp,rsinp,exp1
      real(8) gamma,fax,fay,faz,fcx,fcy,fcz,fb1x,fb1y,fb1z,fd1x,fd1y
      real(8) fd1z,scale,xad,yad,zad,rad,chgprd,coul,fcoul,fx,fy,fz
      real(8) ppp,dlrpot,t1,t2,vk0,vk1,vk2,gk0,gk1,gk2,epsq,engcpe
      real(8) vircpe,rcut,rvdw,engsrp,virsrp,xac,yac,zac,vcon,fcon
      real(8) virc14,b0,rfld0,rfld1,rfld2,alpha,a1,a2,a3,a4,a5,pp,tt
      real(8) cou14_vir,vdw14_vir,strs(6),strs_loc(6)
      real(8), allocatable :: xdab(:),ydab(:),zdab(:)
      real(8), allocatable :: xdbc(:),ydbc(:),zdbc(:)
      real(8), allocatable :: xdcd(:),ydcd(:),zdcd(:)
      data a1,a2,a3/0.254829592d0,-0.284496736d0,1.421413741d0/
      data a4,a5,pp/-1.453152027d0,1.061405429d0,0.3275911d0/
      data fail1,fail2,fail3/0,0,0/

      allocate (xdab(msbad),ydab(msbad),zdab(msbad),stat=fail1)
      allocate (xdbc(msbad),ydbc(msbad),zdbc(msbad),stat=fail2)
      allocate (xdcd(msbad),ydcd(msbad),zdcd(msbad),stat=fail3)
      if(fail1.ne.0.or.fail2.ne.0.or.fail3.ne.0)
     x  call error(idnode,1060)

      twopi=2.d0*pi
      rtwopi=1.d0/twopi
      safe=.true.

c     check size of work arrays

      if((ntdihd-mxnode+1)/mxnode.gt.msbad) call error(idnode,421)

c     block indices

      idih1=(idnode*ntdihd)/mxnode+1
      idih2=((idnode+1)*ntdihd)/mxnode

c     initialise accumulators

      engdih=0.d0
      virdih=0.d0
      dih_fre=0.d0
      dih_vir=0.d0
      strs(:)=0.d0
      strs_loc(:)=0.d0

      if(lsolva)then

        lcomp(3)=.true.
        dih_sol(:)=0.d0
        if(lexcite)dih_exc(:)=0.d0

      endif

      if(keyfce/2.eq.4)then

c     constant terms for shifted coulombic potential

        tt=1.d0/(1.d0+pp*alpha*rcut)
        exp1=exp(-(alpha*rcut)**2)
        vcon=tt*(a1+tt*(a2+tt*(a3+tt*(a4+tt*a5))))*exp1/rcut
        fcon=(vcon+2.d0*(alpha/sqrpi)*exp1)/rcut

      elseif(keyfce/2.eq.5)then

c     constant terms for reaction field potential

        b0=2.d0*(epsq-1.d0)/(2.d0*epsq+1.d0)
        rfld0=b0/rcut**3
        rfld1=(1.d0+b0*0.5d0)/rcut
        rfld2=rfld0*0.5d0
        tt=1.d0/(1.d0+pp*alpha*rcut)
        exp1=exp(-(alpha*rcut)**2)
        vcon=tt*(a1+tt*(a2+tt*(a3+tt*(a4+tt*a5))))*exp1/rcut
        fcon=(vcon+2.d0*(alpha/sqrpi)*exp1)/rcut-rfld0*rcut
        vcon=vcon+rfld2*rcut**2-rfld1

      endif

c     calculate bond vectors

      ii=0
      do i=idih1,idih2

        ii=ii+1

c     indices of bonded atoms

        ia=listdih(ii,2)
        ib=listdih(ii,3)
        ic=listdih(ii,4)
        id=listdih(ii,5)

c     define components of bond vectors

        xdab(ii)=xxx(ia)-xxx(ib)
        ydab(ii)=yyy(ia)-yyy(ib)
        zdab(ii)=zzz(ia)-zzz(ib)

        xdbc(ii)=xxx(ib)-xxx(ic)
        ydbc(ii)=yyy(ib)-yyy(ic)
        zdbc(ii)=zzz(ib)-zzz(ic)

        xdcd(ii)=xxx(ic)-xxx(id)
        ydcd(ii)=yyy(ic)-yyy(id)
        zdcd(ii)=zzz(ic)-zzz(id)

      enddo

c     periodic boundary condition

      call images(imcon,0,1,ii,cell,xdab,ydab,zdab)
      call images(imcon,0,1,ii,cell,xdbc,ydbc,zdbc)
      call images(imcon,0,1,ii,cell,xdcd,ydcd,zdcd)

c     zero dihedral energy accumulator

      engdih=0.d0
      virdih=0.d0

c     zero scaled 1-4 electrostatic and short range potential accumulators

      engc14=0.d0
      virc14=0.d0
      engs14=0.d0
      virs14=0.d0
      vdw14_vir=0.d0
      cou14_vir=0.d0

c     loop over all specified dihedrals

      ii=0
      do i=idih1,idih2

c     define components of bond vectors

        ii=ii+1

        xab=xdab(ii)
        yab=ydab(ii)
        zab=zdab(ii)

        xbc=xdbc(ii)
        ybc=ydbc(ii)
        zbc=zdbc(ii)
        rrbc=1.d0/sqrt(xbc*xbc+ybc*ybc+zbc*zbc)

        xcd=xdcd(ii)
        ycd=ydcd(ii)
        zcd=zdcd(ii)

        xac=xab+xbc
        yac=yab+ybc
        zac=zab+zbc

c     construct first dihedral vector

        pbx=yab*zbc-zab*ybc
        pby=zab*xbc-xab*zbc
        pbz=xab*ybc-yab*xbc
        pb2=pbx*pbx+pby*pby+pbz*pbz
        rpb1=1.d0/sqrt(pb2)
        rpb2=rpb1*rpb1

c     construct second dihedral vector

        pcx=ybc*zcd-zbc*ycd
        pcy=zbc*xcd-xbc*zcd
        pcz=xbc*ycd-ybc*xcd
        pc2=pcx*pcx+pcy*pcy+pcz*pcz
        rpc1=1.d0/sqrt(pc2)
        rpc2=rpc1*rpc1

c     determine dihedral angle

        pbpc=pbx*pcx+pby*pcy+pbz*pcz
        cosp=pbpc*rpb1*rpc1
        sinp=(xbc*(pcy*pbz-pcz*pby)+ybc*(pbx*pcz-pbz*pcx)+
     x    zbc*(pcx*pby-pcy*pbx))*(rpb1*rpc1*rrbc)

        phi=atan2(sinp,cosp)

c     avoid singularity in sinp

        sinp=sign(max(1.d-8,abs(sinp)),sinp)
        rsinp=1.d0/sinp

c     selection of potential energy function type

        kk=listdih(ii,1)

c     calculate potential energy and scalar force term

        if(keydih(kk).eq.1)then

c     key=1 for torsion dihedral potential

          dterm=prmdih(kk,1)*(1.d0+cos(prmdih(kk,3)*phi-
     x      prmdih(kk,2)))
          gamma=-rpb1*rpc1*rsinp*prmdih(kk,1)*prmdih(kk,3)*
     x      sin(prmdih(kk,3)*phi-prmdih(kk,2))

        else if(keydih(kk).eq.2)then

c     key=2 for harmonic improper dihedral

          phi=phi-prmdih(kk,2)
          phi=phi-nint(phi*rtwopi)*twopi
          dterm=0.5d0*prmdih(kk,1)*(phi*phi)
          gamma=rpb1*rpc1*rsinp*prmdih(kk,1)*phi

        else if(keydih(kk).eq.3)then

c     key=3 for harmonic cosine dihedral

          dterm=0.5d0*prmdih(kk,1)*(cos(phi)-
     x      cos(prmdih(kk,2)))**2
          gamma=-rpb1*rpc1*prmdih(kk,1)*(cos(phi)-cos(prmdih(kk,2)))

        else if(keydih(kk).eq.4)then

c     key=4 for 3-term cosine dihedral

          dterm=0.5d0*(prmdih(kk,1)*(1.d0+cos(phi))+
     x      prmdih(kk,2)*(1.d0-cos(2.d0*phi))+prmdih(kk,3)*
     x      (1.d0+cos(3.d0*phi)))
          gamma=-rpb1*rpc1*rsinp*0.5d0*(prmdih(kk,1)*sin(phi)-
     x      2.d0*prmdih(kk,2)*sin(2.d0*phi)+3.d0*prmdih(kk,3)*
     x      sin(3.d0*phi))

        else if(keydih(kk).eq.5)then

c     key=5 for ryckaert-bellemans potential
c     chem.phys.lett., vol.30, p.123, 1975.
c     ATTENTION !!! Modified to have trans configuration
c     correspond to phi=180 rather than
c     phi=0 as in original form.

          dterm=prmdih(kk,1)*(1.116d0-1.462d0*cos(phi)-
     x      1.578d0*(cos(phi))**2+0.368d0*(cos(phi))**3+
     x      3.156d0*(cos(phi))**4+3.788d0*(cos(phi))**5)
          gamma=prmdih(kk,1)*(1.462d0+3.156d0*cos(phi)-
     x      1.104d0*(cos(phi))**2-12.624d0*(cos(phi))**3-
     x      18.94d0*(cos(phi))**4)*rpb1*rpc1

        else if(keydih(kk).eq.6)then

c     key=6 for fluorinated ryckaert-bellemans potential
c     Rice at al., JCP 104, 2101, (1996).

          dterm=prmdih(kk,1)*(3.55d0-2.78d0*cos(phi)-
     x      3.56d0*(cos(phi))**2-1.64d0*(cos(phi))**3+
     x      7.13d0*(cos(phi))**4+12.84d0*(cos(phi))**5+
     x      9.67d0*exp(-56.d0*(phi-pi)**2))
          gamma=(prmdih(kk,1)*(2.78d0+7.12d0*cos(phi)+
     x      4.92d0*(cos(phi))**2-28.52d0*(cos(phi))**3-
     x      64.2d0*(cos(phi))**4)-1083.04d0*(phi-pi)*
     x      exp(-56.0*(phi-pi)**2))*rpb1*rpc1

        else if(keydih(kk).eq.7)then

c     key=7 for opls cosine dihedral

          phi=phi-prmdih(kk,5)
          dterm=prmdih(kk,1)+0.5d0*(prmdih(kk,2)*
     x      (1.d0+cos(phi))+prmdih(kk,3)*(1.d0-cos(2.d0*phi))+
     x      prmdih(kk,4)*(1.d0+cos(3.d0*phi)))
          gamma=-0.5d0*(prmdih(kk,2)*sin(phi)-2.d0*prmdih(kk,3)*
     x      sin(2.d0*phi)+3.d0*prmdih(kk,4)*sin(3.d0*phi))*rpb1*
     x      rpc1*rsinp

        else

c     undefined potential

          safe=.false.
          dterm=0.d0
          gamma=0.d0

        endif

c     indices of bonded atoms

        ia=listdih(ii,2)
        ib=listdih(ii,3)
        ic=listdih(ii,4)
        id=listdih(ii,5)

c     metadynamics local definitions

        if(lmetadyn)then

          idrive=driven(ltype(ia))
          jdrive=driven(ltype(ib))
          kdrive=driven(ltype(ic))
          ldrive=driven(ltype(id))

        endif

c     set selection control for angle potential

        lselect=.true.

        if(lexcite)then

c     selected excitation option

          if((atm_fre(ia).ne.1).and.(atm_fre(ib).ne.1).and.
     x      (atm_fre(ic).ne.1).and.(atm_fre(id).ne.1))then

c     reset selection control

            lselect=(atm_fre(ia)+atm_fre(ib)+atm_fre(ic)+
     x        atm_fre(id).eq.0)

            if(lsolva)then
              dih_exc(atmolt(ia))=dih_exc(atmolt(ia))+dterm
            endif

          endif

        elseif(lfree)then

c     selected free energy option

          if((atm_fre(ia).eq.1).or.(atm_fre(ib).eq.1).or.
     x      (atm_fre(ic).eq.1).or.(atm_fre(id).eq.1))then

c     set hamiltonian mixing parameter

            dih_fre=dih_fre-dterm
            dterm=lambda1*dterm
            gamma=lambda1*gamma

          elseif((atm_fre(ia).eq.2).or.(atm_fre(ib).eq.2).or.
     x        (atm_fre(ic).eq.2).or.(atm_fre(id).eq.2))then

c     set hamiltonian mixing parameter

            dih_fre=dih_fre+dterm
            dterm=lambda2*dterm
            gamma=lambda2*gamma

          endif

        endif

        if(lselect)then

c     calculate potential energy

          engdih=engdih+dterm

c     calculate solvation energy for dihedral term

          if(lsolva)then
            dih_sol(atmolt(ia))=dih_sol(atmolt(ia))+dterm
          endif

c     calculate atomic forces for dihedral term

          fax=gamma*((-pcy*zbc+pcz*ybc)-pbpc*rpb2*(-pby*zbc+pbz*ybc))
          fay=gamma*(( pcx*zbc-pcz*xbc)-pbpc*rpb2*( pbx*zbc-pbz*xbc))
          faz=gamma*((-pcx*ybc+pcy*xbc)-pbpc*rpb2*(-pbx*ybc+pby*xbc))

          fcx=gamma*((-pcy*zab+pcz*yab)-pbpc*rpb2*(-pby*zab+pbz*yab))
          fcy=gamma*(( pcx*zab-pcz*xab)-pbpc*rpb2*( pbx*zab-pbz*xab))
          fcz=gamma*((-pcx*yab+pcy*xab)-pbpc*rpb2*(-pbx*yab+pby*xab))

          fb1x=gamma*((-pby*zcd+pbz*ycd)-pbpc*rpc2*(-pcy*zcd+pcz*ycd))
          fb1y=gamma*(( pbx*zcd-pbz*xcd)-pbpc*rpc2*( pcx*zcd-pcz*xcd))
          fb1z=gamma*((-pbx*ycd+pby*xcd)-pbpc*rpc2*(-pcx*ycd+pcy*xcd))

          fd1x=gamma*((-pby*zbc+pbz*ybc)-pbpc*rpc2*(-pcy*zbc+pcz*ybc))
          fd1y=gamma*(( pbx*zbc-pbz*xbc)-pbpc*rpc2*( pcx*zbc-pcz*xbc))
          fd1z=gamma*((-pbx*ybc+pby*xbc)-pbpc*rpc2*(-pcx*ybc+pcy*xbc))

          fxx(ia)=fxx(ia)+fax
          fyy(ia)=fyy(ia)+fay
          fzz(ia)=fzz(ia)+faz

          fxx(ib)=fxx(ib)-fax-fcx+fb1x
          fyy(ib)=fyy(ib)-fay-fcy+fb1y
          fzz(ib)=fzz(ib)-faz-fcz+fb1z

          fxx(ic)=fxx(ic)+fcx-fb1x-fd1x
          fyy(ic)=fyy(ic)+fcy-fb1y-fd1y
          fzz(ic)=fzz(ic)+fcz-fb1z-fd1z

          fxx(id)=fxx(id)+fd1x
          fyy(id)=fyy(id)+fd1y
          fzz(id)=fzz(id)+fd1z

c     stress tensor for dihedral term

          strs(1)=strs(1)+xab*fax+xbc*(fb1x-fcx)-xcd*fd1x
          strs(2)=strs(2)+yab*fax+ybc*(fb1x-fcx)-ycd*fd1x
          strs(3)=strs(3)+zab*fax+zbc*(fb1x-fcx)-zcd*fd1x
          strs(4)=strs(4)+yab*fay+ybc*(fb1y-fcy)-ycd*fd1y
          strs(5)=strs(5)+yab*faz+ybc*(fb1z-fcz)-ycd*fd1z
          strs(6)=strs(6)+zab*faz+zbc*(fb1z-fcz)-zcd*fd1z

        endif

c     metadynamics local parameters

        if(lmetadyn.and.(idrive.or.jdrive.or.kdrive.or.ldrive))then

c     local energy (no virial)

          eng_loc=eng_loc+dterm

c     local forces

          fxx_loc(ia)=fxx_loc(ia)+fax
          fyy_loc(ia)=fyy_loc(ia)+fay
          fzz_loc(ia)=fzz_loc(ia)+faz

          fxx_loc(ib)=fxx_loc(ib)-fax-fcx+fb1x
          fyy_loc(ib)=fyy_loc(ib)-fay-fcy+fb1y
          fzz_loc(ib)=fzz_loc(ib)-faz-fcz+fb1z

          fxx_loc(ic)=fxx_loc(ic)+fcx-fb1x-fd1x
          fyy_loc(ic)=fyy_loc(ic)+fcy-fb1y-fd1y
          fzz_loc(ic)=fzz_loc(ic)+fcz-fb1z-fd1z

          fxx_loc(id)=fxx_loc(id)+fd1x
          fyy_loc(id)=fyy_loc(id)+fd1y
          fzz_loc(id)=fzz_loc(id)+fd1z

c     local stress tensor

          strs_loc(1)=strs_loc(1)+xab*fax+xbc*(fb1x-fcx)-xcd*fd1x
          strs_loc(2)=strs_loc(2)+yab*fax+ybc*(fb1x-fcx)-ycd*fd1x
          strs_loc(3)=strs_loc(3)+zab*fax+zbc*(fb1x-fcx)-zcd*fd1x
          strs_loc(4)=strs_loc(4)+yab*fay+ybc*(fb1y-fcy)-ycd*fd1y
          strs_loc(5)=strs_loc(5)+yab*faz+ybc*(fb1z-fcz)-ycd*fd1z
          strs_loc(6)=strs_loc(6)+zab*faz+zbc*(fb1z-fcz)-zcd*fd1z

        endif

c     calculate 1-4 dihedral interactions (coulombic and short ranged)
c     assumes 1-4 interactions are in the exclude list

        kk=listdih(ii,1)

c     bypass OPLS 1-4 terms (not present)

        if(keydih(kk).ne.7)then

c     1-4 electrostatics : adjust by weighting factor

          scale=prmdih(kk,4)

          xad=xac+xcd
          yad=yac+ycd
          zad=zac+zcd

          rad=sqrt(xad**2+yad**2+zad**2)

c     scaled charge product*dielectric

          chgprd=scale*chge(ia)*chge(id)*r4pie0
          coul=0.d0
          fcoul=0.d0

c     truncation of potential for all schemes except ewald sum

          if(abs(chgprd).gt.1.d-10.and.keyfce.gt.0)then

c     electrostatics by ewald sum

            if(keyfce/2.eq.1.or.keyfce/2.eq.6.or.keyfce/2.eq.7)then

              coul=chgprd/(epsq*rad)
              fcoul=coul/(rad**2)

c     distance dependent dielectric

            elseif(rcut.gt.rad)then

              if(keyfce/2.eq.2)then

                coul=chgprd/(epsq*rad**2)
                fcoul=2.0d0*coul/(rad**2)

c     unmodified coulombic

              else if(keyfce/2.eq.3)then

                coul=chgprd/(epsq*rad)
                fcoul=coul/(rad**2)

c     truncated and shifted coulombic

              else if(keyfce/2.eq.4)then

                tt=1.d0/(1.d0+pp*alpha*rad)
                exp1=exp(-(alpha*rad)**2)
                erc=tt*(a1+tt*(a2+tt*(a3+tt*(a4+tt*a5))))*exp1/rad
                fer=(erc+2.d0*(alpha/sqrpi)*exp1)/rad**2

c     calculate potential energy and forces

                coul=chgprd*(erc-vcon+fcon*(rad-rcut))/epsq
                fcoul=chgprd*(fer-fcon/rad)/epsq

c     reaction field

              else if(keyfce/2.eq.5)then

                tt=1.d0/(1.d0+pp*alpha*rad)
                exp1=exp(-(alpha*rad)**2)
                erc=tt*(a1+tt*(a2+tt*(a3+tt*(a4+tt*a5))))*exp1/rad
                fer=(erc+2.d0*(alpha/sqrpi)*exp1)/rad**2
                coul=chgprd*(erc-vcon+fcon*(rad-rcut)+
     x            rfld2*rad*rad-rfld1)
                fcoul=chgprd*(fer-fcon/rad-rfld0)

              elseif(keyfce/2.eq.0)then

                coul=0.d0
                fcoul=0.d0

              else

                call error(idnode,446)

              endif

            endif

c     set selection control for coulombic 1-4 terms

            lselect=.true.

c     set double index

            if(lsolva)kkk=loc2(atmolt(ia),atmolt(id))

            if(lexcite)then

c     selected excitation option

              if((atm_fre(ia).ne.1).and.(atm_fre(id).ne.1))then

c     set selection control

                lselect=(atm_fre(ia)+atm_fre(id).eq.0)

                if(lsolva)cou_exc(kkk)=cou_exc(kkk)+coul

              endif

            elseif(lfree)then

c     selected free energy option

              if((atm_fre(ia).eq.1).or.(atm_fre(id).eq.1))then

c     set hamiltonian mixing parameter

                cou_fre=cou_fre-coul
                cou14_vir=cou14_vir+fcoul*rad**2
                coul=lambda1*coul
                fcoul=lambda1*fcoul

              elseif((atm_fre(ia).eq.2).or.(atm_fre(id).eq.2))then

c     set hamiltonian mixing parameter

                cou_fre=cou_fre+coul
                cou14_vir=cou14_vir-fcoul*rad**2
                coul=lambda2*coul
                fcoul=lambda2*fcoul

              endif

            endif

            if(lselect)then

c     electrostatic energy and virial for 1-4 term

              engc14=engc14+coul
              virc14=virc14-fcoul*rad**2

c     solvation energy for coulombic 1-4 term

              if(lsolva)cou_sol(kkk)=cou_sol(kkk)+coul

c     coulombic force for 1-4 term

              fx=fcoul*xad
              fy=fcoul*yad
              fz=fcoul*zad

              fxx(ia)=fxx(ia)+fx
              fyy(ia)=fyy(ia)+fy
              fzz(ia)=fzz(ia)+fz

              fxx(id)=fxx(id)-fx
              fyy(id)=fyy(id)-fy
              fzz(id)=fzz(id)-fz

c     stress tensor for coulombic 1-4 term

              strs(1)=strs(1)+xad*fx
              strs(2)=strs(2)+xad*fy
              strs(3)=strs(3)+xad*fz
              strs(4)=strs(4)+yad*fy
              strs(5)=strs(5)+yad*fz
              strs(6)=strs(6)+zad*fz

            endif

c     metadynamics local parameters

            if(lmetadyn.and.(idrive.or.ldrive))then

c     local energy and virial

              eng_loc=eng_loc+coul
              vir_loc=vir_loc-fcoul*rad**2

c     local forces

              fxx_loc(ia)=fxx_loc(ia)+fx
              fyy_loc(ia)=fyy_loc(ia)+fy
              fzz_loc(ia)=fzz_loc(ia)+fz

              fxx_loc(id)=fxx_loc(id)-fx
              fyy_loc(id)=fyy_loc(id)-fy
              fzz_loc(id)=fzz_loc(id)-fz

c     local stress tensor

              strs_loc(1)=strs_loc(1)+xad*fx
              strs_loc(2)=strs_loc(2)+xad*fy
              strs_loc(3)=strs_loc(3)+xad*fz
              strs_loc(4)=strs_loc(4)+yad*fy
              strs_loc(5)=strs_loc(5)+yad*fz
              strs_loc(6)=strs_loc(6)+zad*fz

            endif

          endif

c     1-4 short ranged : adjust by weighting factor

          scale=prmdih(kk,5)
          gamma=0.d0

          if(mod(keyfce,2).eq.1)then

c     atomic and potential function indices

            ka=max(ltype(ia),ltype(id))
            kb=min(ltype(ia),ltype(id))
            k=lstvdw((ka*(ka-1))/2+kb)

            if(abs(scale*vvv(1,k)).gt.1.d-10)then

c     apply truncation of potential

              if(rvdw.gt.rad)then

c     determine interpolation panel for force arrays

                l=int(rad/dlrpot)
                ppp=rad/dlrpot-dble(l)

c     calculate interaction energy using 3-point interpolation

                vk0=vvv(l,k)
                vk1=vvv(l+1,k)
                vk2=vvv(l+2,k)

                t1=vk0+(vk1-vk0)*ppp
                t2=vk1+(vk2-vk1)*(ppp-1.0d0)

                srpot=scale*(t1+(t2-t1)*ppp*0.5d0)

c     calculate forces using 3-point interpolation

                gk0=ggg(l,k)
                gk1=ggg(l+1,k)
                gk2=ggg(l+2,k)

                t1=gk0+(gk1-gk0)*ppp
                t2=gk1+(gk2-gk1)*(ppp-1.0d0)

                gamma=scale*(t1+(t2-t1)*ppp*0.5d0)/(rad**2)

c     set selection control for short ranged 1-4 terms

                lselect=.true.

c     set double index

                if(lsolva)kkk=loc2(atmolt(ia),atmolt(id))

                if(lexcite)then

c     selected excitation option

                  if((atm_fre(ia).ne.1).and.(atm_fre(id).ne.1))then

c     set selection control

                    lselect=(atm_fre(ia)+atm_fre(id).eq.0)

                    if(lsolva)vdw_exc(kkk)=vdw_exc(kkk)+srpot

                  endif

                elseif(lfree)then

c     selected free energy option

                  if((atm_fre(ia).eq.1).or.(atm_fre(id).eq.1))then

c     set hamiltonian mixing parameter

                    vdw_fre=vdw_fre-srpot
                    vdw14_vir=vdw14_vir+gamma*rad**2
                    srpot=lambda1*srpot
                    gamma=lambda1*gamma

                  elseif((atm_fre(ia).eq.2).or.(atm_fre(id).eq.2))then

c     set hamiltonian mixing parameter

                    vdw_fre=vdw_fre+srpot
                    vdw14_vir=vdw14_vir-gamma*rad**2
                    srpot=lambda2*srpot
                    gamma=lambda2*gamma

                  endif

                endif

                if(lselect)then

c     short ranged energy and virial for 1-4 term

                  engs14=engs14+srpot
                  virs14=virs14-gamma*rad**2

c     solvation energy for short ranged 1-4 term

                  if(lsolva)vdw_sol(kkk)=vdw_sol(kkk)+srpot

c     short ranged forces for 1-4 term

                  fx=gamma*xad
                  fy=gamma*yad
                  fz=gamma*zad

                  fxx(ia)=fxx(ia)+fx
                  fyy(ia)=fyy(ia)+fy
                  fzz(ia)=fzz(ia)+fz

                  fxx(id)=fxx(id)-fx
                  fyy(id)=fyy(id)-fy
                  fzz(id)=fzz(id)-fz

c     stress tensor for short ranged 1-4 term

                  strs(1)=strs(1)+xad*fx
                  strs(2)=strs(2)+xad*fy
                  strs(3)=strs(3)+xad*fz
                  strs(4)=strs(4)+yad*fy
                  strs(5)=strs(5)+yad*fz
                  strs(6)=strs(6)+zad*fz

                endif

c     metadynamics local parameters

                if(lmetadyn.and.(idrive.or.ldrive))then

c     local energy and virial

                  eng_loc=eng_loc+srpot
                  vir_loc=vir_loc-gamma*rad**2

c     local forces

                  fxx_loc(ia)=fxx_loc(ia)+fx
                  fyy_loc(ia)=fyy_loc(ia)+fy
                  fzz_loc(ia)=fzz_loc(ia)+fz

                  fxx_loc(id)=fxx_loc(id)-fx
                  fyy_loc(id)=fyy_loc(id)-fy
                  fzz_loc(id)=fzz_loc(id)-fz

c     stress tensor for short ranged 1-4 term

                  strs_loc(1)=strs_loc(1)+xad*fx
                  strs_loc(2)=strs_loc(2)+xad*fy
                  strs_loc(3)=strs_loc(3)+xad*fz
                  strs_loc(4)=strs_loc(4)+yad*fy
                  strs_loc(5)=strs_loc(5)+yad*fz
                  strs_loc(6)=strs_loc(6)+zad*fz

                endif

              endif

            endif

          endif

        endif

      enddo

c     complete stress tensor

      stress(1)=stress(1)+strs(1)
      stress(2)=stress(2)+strs(2)
      stress(3)=stress(3)+strs(3)
      stress(4)=stress(4)+strs(2)
      stress(5)=stress(5)+strs(4)
      stress(6)=stress(6)+strs(5)
      stress(7)=stress(7)+strs(3)
      stress(8)=stress(8)+strs(5)
      stress(9)=stress(9)+strs(6)

      if(lmetadyn)then

        stress_loc(1)=stress_loc(1)+strs_loc(1)
        stress_loc(2)=stress_loc(2)+strs_loc(2)
        stress_loc(3)=stress_loc(3)+strs_loc(3)
        stress_loc(4)=stress_loc(4)+strs_loc(2)
        stress_loc(5)=stress_loc(5)+strs_loc(4)
        stress_loc(6)=stress_loc(6)+strs_loc(5)
        stress_loc(7)=stress_loc(7)+strs_loc(3)
        stress_loc(8)=stress_loc(8)+strs_loc(5)
        stress_loc(9)=stress_loc(9)+strs_loc(6)

      endif

c     sum contributions to potentials

      if(mxnode.gt.1)then

        buffer(1)=engdih
        buffer(2)=engc14
        buffer(3)=virc14
        buffer(4)=engs14
        buffer(5)=virs14
        buffer(6)=dih_fre
        buffer(7)=dih_vir
        buffer(8)=vdw14_vir
        buffer(9)=cou14_vir
        call gdsum(buffer(1),9,buffer(10))
        engdih=buffer(1)
        engc14=buffer(2)
        virc14=buffer(3)
        engs14=buffer(4)
        virs14=buffer(5)
        dih_fre=buffer(6)
        dih_vir=buffer(7)
        vdw14_vir=buffer(8)
        cou14_vir=buffer(9)

        if(lsolva)then

          call gdsum(dih_sol(1),mxtmls,buffer(1))
          if(lexcite)call gdsum(dih_exc(1),mxtmls,buffer(1))

        endif

      endif

      engcpe=engcpe+engc14
      vircpe=vircpe+virc14
      engsrp=engsrp+engs14
      virsrp=virsrp+virs14
      vdw_vir=vdw_vir+vdw14_vir
      cou_vir=cou_vir+cou14_vir

c     check for undefined potentials

      if(mxnode.gt.1)call gstate(safe)
      if(.not.safe)call error(idnode,448)

      deallocate (xdab,ydab,zdab,stat=fail1)
      deallocate (xdbc,ydbc,zdbc,stat=fail2)
      deallocate (xdcd,ydcd,zdcd,stat=fail3)

      return
      end subroutine dihfrc

      end module dihedral_module


text:
      subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
c  subroutine splder evaluates in a number of points x(i),i=1,2,...,m
c  the derivative of order nu of a spline s(x) of degree k,given in
c  its b-spline representation.
c
c  calling sequence:
c     call splder(t,n,c,k,nu,x,y,m,wrk,ier)
c
c  input parameters:
c    t    : array,length n, which contains the position of the knots.
c    n    : integer, giving the total number of knots of s(x).
c    c    : array,length n, which contains the b-spline coefficients.
c    k    : integer, giving the degree of s(x).
c    nu   : integer, specifying the order of the derivative. 0<=nu<=k
c    x    : array,length m, which contains the points where the deriv-
c           ative of s(x) must be evaluated.
c    m    : integer, giving the number of points where the derivative
c           of s(x) must be evaluated
c    wrk  : real array of dimension n. used as working space.
c
c  output parameters:
c    y    : array,length m, giving the value of the derivative of s(x)
c           at the different points.
c    ier  : error flag
c      ier = 0 : normal return
c      ier =10 : invalid input data (see restrictions)
c
c  restrictions:
c    0 <= nu <= k
c    m >= 1
c    t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
c
c  other subroutines required: fpbspl
c
c  references :
c    de boor c : on calculating with b-splines, j. approximation theory
c                6 (1972) 50-62.
c    cox m.g.  : the numerical evaluation of b-splines, j. inst. maths
c                applics 10 (1972) 134-149.
c   dierckx p. : curve and surface fitting with splines, monographs on
c                numerical analysis, oxford university press, 1993.
c
c  author :
c    p.dierckx
c    dept. computer science, k.u.leuven
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c  latest update : march 1987
c
c  ..scalar arguments..
      integer n,k,nu,m,ier
c  ..array arguments..
      real t(n),c(n),x(m),y(m),wrk(n)
c  ..local scalars..
      integer i,j,kk,k1,k2,l,ll,l1,l2,nk1,nk2,nn
      real ak,arg,fac,sp,tb,te
c  ..local arrays ..
      real h(6)
c  before starting computations a data check is made. if the input data
c  are invalid control is immediately repassed to the calling program.
      ier = 10
      if(nu.lt.0 .or. nu.gt.k) go to 200
      if(m-1) 200,30,10
  10  do 20 i=2,m
        if(x(i).lt.x(i-1)) go to 200
  20  continue
  30  ier = 0
c  fetch tb and te, the boundaries of the approximation interval.
      k1 = k+1
      nk1 = n-k1
      tb = t(k1)
      te = t(nk1+1)
c  the derivative of order nu of a spline of degree k is a spline of
c  degree k-nu,the b-spline coefficients wrk(i) of which can be found
c  using the recurrence scheme of de boor.
      l = 1
      kk = k
      nn = n
      do 40 i=1,nk1
         wrk(i) = c(i)
  40  continue
      if(nu.eq.0) go to 100
      nk2 = nk1
      do 60 j=1,nu
         ak = kk
         nk2 = nk2-1
         l1 = l
         do 50 i=1,nk2
            l1 = l1+1
            l2 = l1+kk
            fac = t(l2)-t(l1)
            if(fac.le.0.) go to 50
            wrk(i) = ak*(wrk(i+1)-wrk(i))/fac
  50     continue
         l = l+1
         kk = kk-1
  60  continue
      if(kk.ne.0) go to 100
c  if nu=k the derivative is a piecewise constant function
      j = 1
      do 90 i=1,m
         arg = x(i)
  70     if(arg.lt.t(l+1) .or. l.eq.nk1) go to 80
         l = l+1
         j = j+1
         go to 70
  80     y(i) = wrk(j)
  90  continue
      go to 200
 100  l = k1
      l1 = l+1
      k2 = k1-nu
c  main loop for the different points.
      do 180 i=1,m
c  fetch a new x-value arg.
        arg = x(i)
        if(arg.lt.tb) arg = tb
        if(arg.gt.te) arg = te
c  search for knot interval t(l) <= arg < t(l+1)
 140    if(arg.lt.t(l1) .or. l.eq.nk1) go to 150
        l = l1
        l1 = l+1
        go to 140
c  evaluate the non-zero b-splines of degree k-nu at arg.
 150    call fpbspl(t,n,kk,arg,l,h)
c  find the value of the derivative at x=arg.
        sp = 0.
        ll = l-k1
        do 160 j=1,k2
          ll = ll+1
          sp = sp+wrk(ll)*h(j)
 160    continue
        y(i) = sp
 180  continue
 200  return
      end

ext: f

max_stars_repo_name: harbik/dierckx-sys

max_issues_repo_name: harbik/dierckx-sys

max_forks_repo_name: harbik/dierckx-sys


text:
C
C
*+ plot_cursor_read

       subroutine plot_cursor_read(map_array,report,s)
C      -----------------------------------------------
C
C Read cursor position and return UV, RA/DEC and value
C
C Given:
C  map array
       real*4    map_array(*)
C  flag requesting results are summarized on output
       logical   report
C
C Updated:
C  error status
       integer   s
C
C This is a high level routine and may be called at any time. The routine
C checks that a valid plot is open. The output consists of:
C
C    the nearest UV pixel and value %u, %v, %uv-pos, %uv-value
C    the interpolated position and value %map-u, %map-v, %map-pos, %map-value
C    primary-beam corrected values %map-pbvalue
C    RA and DEC at cursor location %ra, %dec
C
C These data are returned via standard named parameters and if report is
C set to true then they are also reported on the output unit.
C
*-
       include '../include/plt_basic_defn.inc'
       include '../include/plt_image_defn.inc'
       include '/mrao/include/constants.inc'
       include '/mrao/include/chrlib_functions.inc'

C graphic position
       real*4       gx, gy
       character    ch_typed*1, rnumber*20, value*80
       integer      i1, i2, len_v

C cursor position in map coordinates
       real*4       pos_u, pos_v
       real*8       dpos_u, dpos_v

C value at nearest pixel, interpolated value, and primary-beam corrected
       real*4       value_int, value_pix, pb_value
       character    string_int*14, string_pix*14, string_pb*14
       integer      len_val_int, len_val_pix, len_val_pb

C RA/DEC of current position
       real*8       ra, dec, ra0, dec0, epoch, rapnt, decpnt
       real*4       freq
       integer      tscope

C parameters to hold RA/DEC position
       integer      i_hour, i_min_time, i_deg, i_min_arc
       real*4       sec_time, sec_arc

C functions and flags
       integer      pgcurse, iout
       logical      cursor_available, cmd_dblev

       if (s.ne.0) return

       if (.not.plot_open) then
         call cmd_wrerr('CURSOR-POSITION',
     *                  'Plot not open - unable to access cursor')
         return
       end if
       if (.not.frame_init) then
         call cmd_wrerr('CURSOR-POSITION',
     *                  'Frame not initialised for image')
         return
       end if
       call pgqcur(cursor_available)
       if (.not.cursor_available) then
         call cmd_wrerr('CURSOR-POSITION',
     *                   'Device does not support input')
         return
       end if

C find initial position for cursor
       if (cmd_dblev(1)) then
         call io_wrout('.. move cursor: type any character to end')
       endif

       pos_u = (frame_win(1)+frame_win(2))/2.0
       pos_v = (frame_win(3)+frame_win(4))/2.0

C read cursor position
       s = pgcurse(pos_u,pos_v,ch_typed) - 1
       if (s.ne.0) then
         call cmd_wrerr('CURSOR-POSITION','Error reading CURSOR')
         return
       end if

C set value for graphic position and character
       gx = (pos_u-frame_win(1))*(frame_vp(2)-frame_vp(1))/
     *      (frame_win(2)-frame_win(1)) + frame_vp(1)
       gy = (pos_v-frame_win(3))*(frame_vp(4)-frame_vp(3))/
     *      (frame_win(4)-frame_win(3)) + frame_vp(3)
       call cmd_setlocal( '%gr-char', ch_typed, s )
       call chr_chrtoc(gx,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal( '%gr-x', rnumber(i1:i2), s )
       call chr_chrtoc(gy,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal( '%gr-y', rnumber(i1:i2), s )

C determine map value at cursor position and validity
       dpos_u = pos_u
       dpos_v = pos_v
       call ruvval(map_array,dpos_u,dpos_v,2,value_int,s)
       s = 0
       call iuvval(map_array,nint(pos_u),nint(pos_v),value_pix,s)
       call chr_chrtoc(value_int,string_int,len_val_int)
       call chr_chrtoc(value_pix,string_pix,len_val_pix)
       call cmd_setlocal('%map-value',string_int(1:len_val_pix),s)
       call cmd_setlocal('%uv-value',string_pix(1:len_val_pix),s)
       call chr_chitoc(nint(pos_u),rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%u',rnumber(i1:i2),s)
       value = rnumber(i1:i2)
       call chr_chitoc(nint(pos_v),rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%v',rnumber(i1:i2),s)
       len_v = chr_lenb(value)
       value = value(1:len_v)//','//rnumber(i1:i2)
       len_v = chr_lenb(value)
       call cmd_setlocal( '%uv-pos', value(1:len_v), s )

       call chr_chrtoc(pos_u,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%map-u',rnumber(i1:i2),s)
       value = rnumber(i1:i2)
       call chr_chrtoc(pos_v,rnumber,i2)
       i1 = chr_intlc(rnumber)
       call cmd_setlocal('%map-v',rnumber(i1:i2),s)
       len_v = chr_lenb(value)
       value = value(1:len_v)//','//rnumber(i1:i2)
       len_v = chr_lenb(value)
       call cmd_setlocal( '%map-pos', value(1:len_v), s )

C output results if reporting
       if (report) then
         call io_enqout(iout)
         write(iout,10) nint(pos_u),nint(pos_v),value_pix
         write(iout,11) pos_u,pos_v,value_int
10       format(' Map pixel    : ',2I6,12X,'  Pixel-value : ',1PE12.3)
11       format(' Map position : ',2F12.3, '  Map-value   : ',1PE12.3)
       endif

C return if point not on map
       if (s.ne.0) then
         s = 0
         return
       endif

C find primary beam corrected value
       call entpnt(tscope, freq, rapnt, decpnt, s)
       call enepoch( epoch, s )
       call uvtord(dpos_u,dpos_v,ra0,dec0,s)
       call precrd(epoch,ra0,dec0,1950.0d0,ra,dec)
       call pbcorr(ra0,dec0,rapnt,decpnt,tscope,pb_value,s)
       if (pb_value.gt.0.0) then
         value_int = value_int/pb_value
         call chr_chrtoc(value_int,string_pb,len_val_pb)
       else
         call chr_chrtoc(0.0,string_pb,len_val_pb)
       end if
       call cmd_setlocal('%map-pbvalue',string_pb(1:len_val_pb),s)
       if (report) then
         write(iout,12)value_int
12       format('                ',24X, '  PB-value    : ',1PE12.3)

       endif

C find RA/DEC
       call frrads('HMS',ra,i_hour,i_min_time,sec_time)
       call frrads('DMS',dec,i_deg,i_min_arc,sec_arc)
       value = ' '
       write(value,1) i_hour,i_min_time,sec_time
       len_v = chr_lenb(value)
       call cmd_setlocal( '%RA', value(1:len_v), s )
 1     format(I3,1X,I2,1X,F7.3)
       value = ' '
       write(value,2) i_deg,i_min_arc,sec_arc
       len_v = chr_lenb(value)
       call cmd_setlocal( '%DEC', value(1:len_v), s )
 2     format(I3,1X,I2,1X,F7.2)

       if (report) then
         write(iout,13) i_hour,i_min_time,sec_time,
     *                  i_deg,i_min_arc,sec_arc
13       format(' RA DEC       : ',I4,1X,I2,1X,F7.3,I4,1X,I2,1X,F7.3)
       endif

       call cmd_err( s,'PLOT_CURSOR_READ',' ')

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
      subroutine wricon2(lundia,lun   ,bndval,itypc ,nobnd ,notims,
     *                   kmax  ,lstci ,nocon ,
     *                   tstart,dtmin ,
     *                   namcon,nambnd,itdate              )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: wricon2.f 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/tools_gpl/nesthd2/packages/nesthd2/src/wricon2.f $
!***********************************************************************
! Deltares                         marine and coastal management
!
! subroutine         : wricon2
! version            : v1.0
! date               : September 1997
! programmer         : Theo van der Kaaij
!
! function           : Writes transport bc. to TRISULA attribute file
!                      (new format)
! limitations        :
! subroutines called :
!***********************************************************************

      character(len=37), dimension(13) :: fmtbcc
      character(len=20), dimension(16) :: keywrd

      integer      itdate

      real         bndval (nobnd ,notims,kmax  ,lstci ,2)

      character*1  eol
      character*1  quote
      character*3  cdum
      character*17 layer
      character*20 nambnd (nobnd )
      character*20 namcon (lstci )
      character*40 cntain

      data         fmtbcc/'(a1,i5,                       t89,a1)',
     *                    '(a20,a1,a63,a3,a1,            t89,a1)',
     *                    '(a20,a1,a10,a1,a40,           t89,a1)',
     *                    '(2(a20,a1),                   t89,a1)',
     *                    '(a20,                         t89,a1)',
     *                    '(a20,                         t89,a1)',
     *                    '(a20,a1,a15,a1,               t89,a1)',
     *                    '(a20,i9,                      t89,a1)',
     *                    '(a20,a1,a7,a1,                t89,a1)',
     *                    '(a20,a1,a6,a1,                t89,a1)',
     *                    '(a20,a1,a37,a1,1x,2(a10,a1),  t89,a1)',
     *                    '(a20,i6,                      t89,a1)',
     *                    '(f12.4,    5ES14.6,           t89,a1)'/

      data         keywrd/'table-name          ',
     *                    'contents            ',
     *                    'location            ',
     *                    'geo-coordinates     ',
     *                    'metric-coordinates  ',
     *                    'frequency           ',
     *                    'frequency-unit      ',
     *                    'time-function       ',
     *                    'reference-time      ',
     *                    'time-unit           ',
     *                    'time-step           ',
     *                    'interpolation       ',
     *                    'extrapolation       ',
     *                    'parameter           ',
     *                    'unit                ',
     *                    'records-in-table    '/
!----------------------------------------------------------------------
!     Initialisation
!----------------------------------------------------------------------
      write (*     ,
     *'('' >>> Writing transport boundary conditions <<<'')')
      write (lundia,
     *'('' >>> Writing transport boundary conditions <<<'')')

      cdum = '   '
      quote  = CHAR (39)
      eol    = ' '
!     eol    = CHAR (10)
!
      cntain = ' # at ends A&B of open boundary segment '

      do 10 ifmt = 1, 13
         write (fmtbcc(ifmt)(29:29),'(a1)') 't'
         write (fmtbcc(ifmt)(30:33),'(i4.4)') max(12 + 2*kmax*14 + 1,89)
   10 continue

!----------------------------------------------------------------------
!     Write dummy record to avoid binary characters on first record (NT)
!----------------------------------------------------------------------

      write (lun,fmtbcc(5)) ' ',eol
      backspace (lun)
!
!----------------------------------------------------------------------
!     Uniform and linear boundary conditions
!     (uniform also a linear profile with identical surface and bed
!      values)
!----------------------------------------------------------------------
      do 100 ibnd = 1, nobnd
         do 100 ic = 1, nocon
            do 100 itim = 1, notims
!               irec   = (ibnd-1)*(lstci*(notims+11)) +
!    *                            (ic-1)*(notims+11)  + 1
!---------------------------------------------------------------------
!---------------Write first 11 records with block information for
!               time IB=1
!               PARNAM starts at 2 (first id for DIS)
!---------------------------------------------------------------------
                if (itim .eq. 1) then
                  write(lun   ,fmtbcc( 2)) keywrd( 1),
     *                  quote,'Nested transport bc '//namcon(ic)//
     *                  ' '//nambnd(ibnd)//'  ',
     *                  cdum,quote,eol
                  if (itypc .eq. 1) then
                     write(lun   ,fmtbcc( 3)) keywrd( 2),
     *                     quote,'Uniform   '    ,quote,cntain,eol
                  elseif (itypc .eq. 2) then
                     write(lun   ,fmtbcc( 3)) keywrd( 2),
     *                     quote,'Linear    '    ,quote,cntain,eol
                  else
                     write(lun   ,fmtbcc( 3)) keywrd( 2),
     *                     quote,'3d-profile',quote,cntain,eol
                  endif
                  write(lun   ,fmtbcc( 4)) keywrd( 3),
     *                  quote,nambnd(ibnd),quote,eol
!                 write(lun   ,fmtbcc( 5)) keywrd( 4),eol
!                 write(lun   ,fmtbcc( 6)) keywrd( 5),eol
                  write(lun   ,fmtbcc( 7)) keywrd( 8),
     *                  quote,'non-equidistant',quote,eol
                  write(lun   ,fmtbcc( 8)) keywrd( 9),
     *                  itdate,eol
                  write(lun   ,fmtbcc( 9)) keywrd(10),
     *                  quote,'minutes',quote,eol
                  write(lun   ,fmtbcc(10)) keywrd(12),
     *                  quote,'linear',quote,eol
                  if (itypc .ne. 3) then
                     write(lun   ,fmtbcc(11)) keywrd(14),
     *                  quote,'time                '//
     *                        '                 ',
     *                  quote,keywrd(15),
     *                  quote,' ',quote,eol
                     write(lun   ,fmtbcc(11)) keywrd(14),
     *                  quote,namcon(ic)//' End A surface   ',
     *                  quote,keywrd(15),
     *                  quote,' ',quote,eol
                     if (itypc .eq. 2) then
                        write(lun   ,fmtbcc(11)) keywrd(14),
     *                     quote,namcon(ic)//' End A bed       ',
     *                     quote,keywrd(15),
     *                     quote,' ',quote,eol
                     endif
                     write(lun   ,fmtbcc(11)) keywrd(14),
     *                  quote,namcon(ic)//' End B surface   ',
     *                  quote,keywrd(15),
     *                  quote,' ',quote,eol
                     if (itypc .eq. 2) then
                        write(lun   ,fmtbcc(11)) keywrd(14),
     *                     quote,namcon(ic)//' End B bed       ',
     *                     quote,keywrd(15),
     *                     quote,' ',quote,eol
                     endif
                  else
                     write(lun   ,fmtbcc(11)) keywrd(14),
     *                  quote,'time                '//
     *                        '                 ',
     *                  quote,keywrd(15),
     *                  quote,' ',quote,eol
                     do 110 k = 1, kmax
                        layer = ' End A layer:    '
                        write (layer(15:17),'(a3)') cdum
                        write (layer(15:17),'(i3)') k
                        write(lun   ,fmtbcc(11)) keywrd(14),
     *                     quote,namcon(ic)//layer,
     *                     quote,keywrd(15),
     *                     quote,' ',quote,eol
  110                continue
                     do 120 k = 1, kmax
                        layer = ' End B layer:    '
                        write (layer(15:17),'(a3)') cdum
                        write (layer(15:17),'(i3)') k
                        write(lun   ,fmtbcc(11)) keywrd(14),
     *                     quote,namcon(ic)//layer,
     *                     quote,keywrd(15),
     *                     quote,' ',quote,eol
  120                continue
                  endif
                  write(lun   ,fmtbcc(12)) keywrd(16),
     *                  notims,eol
                endif
!----------------------------------------------------------------------
!---------------Write time depended data to block for constituent
!               L skipping first 11 records with block info
!----------------------------------------------------------------------
!               irec    = irec   + itim   + 11
                if     (itypc .eq. 1) then
                  write (fmtbcc(13) (10:12),'(i3.3)') 2
                  write (lun   ,fmtbcc(13))
     *                   tstart + (itim - 1)*dtmin,
     *                   bndval(ibnd,itim,1   ,ic,1),
     *                   bndval(ibnd,itim,1   ,ic,2),eol
                elseif     (itypc .eq. 2) then
                  write (fmtbcc(13) (10:12),'(i3.3)') 4
                  write (lun   ,fmtbcc(13))
     *                   tstart + (itim - 1)*dtmin,
     *                   bndval(ibnd,itim,1   ,ic,1),
     *                   bndval(ibnd,itim,kmax,ic,1),
     *                   bndval(ibnd,itim,1   ,ic,2),
     *                   bndval(ibnd,itim,kmax,ic,2),eol
                else
                   write (fmtbcc(13) (10:12),'(i3.3)') 2*kmax
                   write (lun   ,fmtbcc(13))
     *                   tstart + (itim - 1)*dtmin,
     *                  (bndval(ibnd,itim,k,ic,1),k=1,kmax),
     *                  (bndval(ibnd,itim,k,ic,2),k=1,kmax),eol
                endif
  100 continue

      return

      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
! { dg-do compile }
! { dg-options "-O3 -ffast-math" }

      SUBROUTINE ZUNG2L( M, N, K, A, LDA, TAU, WORK, INFO )
      COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
      IF( M.LT.0 ) THEN
      END IF
         CALL ZLARF( 'LEFT', M-N+II, II-1, A( 1, II ), 1, TAU( I ), A,
     $               LDA, WORK )
         CALL ZSCAL( M-N+II-1, -TAU( I ), A( 1, II ), 1 )
         A( M-N+II, II ) = ONE - TAU( I )
      END
      SUBROUTINE ZLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
      CHARACTER         SIDE(*)
      LOGICAL            LSAME
      COMPLEX*16         C( LDC, * ), V(*), WORK(*), TAU
      IF( LSAME( SIDE, 'L' ) ) THEN
         IF( TAU.NE.ZERO ) THEN
            CALL ZGEMV( 'CONJUGATE TRANSPOSE', M, N, ONE, C, LDC, V,
     $                  INCV, ZERO, WORK, 1 )
         END IF
      END IF
      END
      LOGICAL          FUNCTION LSAME( CA, CB )
      CHARACTER         CA(*), CB(*)
      END
      SUBROUTINE ZGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
     $                   BETA, Y, INCY )
      COMPLEX*16         A( LDA, * ), X( * ), Y( * )
      CHARACTER         TRANS(*)
      LOGICAL            LSAME
      IF( LSAME( TRANS, 'N' ) )THEN
         IF( INCY.EQ.1 )THEN
               IF( X( JX ).NE.ZERO )THEN
                     Y( I ) = Y( I ) + TEMP*A( I, J )
               END IF
         END IF
      END IF
      END
      SUBROUTINE  ZSCAL(N,ZA,ZX,INCX)
      COMPLEX*16 ZA,ZX(1)
      IF( N.LE.0 .OR. INCX.LE.0 )RETURN
   20 DO 30 I = 1,N
        ZX(I) = ZA*ZX(I)
   30 CONTINUE
      END


  !
  !  Copyright (c) 2006-2015, The Regents of the University of California,
  !  through Lawrence Berkeley National Laboratory (subject to receipt of any
  !  required approvals from the U.S. Dept. of Energy) and the Paul Scherrer
  !  Institut (Switzerland).  All rights reserved.!
  !
  !  License: see file COPYING in top level of source distribution.
  !
  include 'H5hut.f90'

  program write_step_attribs
    use H5hut
    implicit none

#if defined(PARALLEL_IO)
    include 'mpif.h'
#endif

    integer*8, parameter :: verbosity_level =          1
    character (len=*), parameter :: FNAME =            "example_step_attribs.h5"

    character (len=*), parameter :: ATTR_STRING =      "StepAttrString"
    character (len=*), parameter :: ATTR_I4 =          "StepAttrInt32"
    character (len=*), parameter :: ATTR_I8 =          "StepAttrInt64"
    character (len=*), parameter :: ATTR_R4 =          "StepAttrFloat32"
    character (len=*), parameter :: ATTR_R8 =          "StepAttrFloat64"

    character (len=*),parameter :: string_value =      "This is a string attribute bound to this step."

    integer*4, parameter, dimension(*) :: i4_value =   (/0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144/)
    integer*8, parameter, dimension(*) :: i8_value =   (/42, 43, 44, 45/)
    real*4, parameter, dimension(*) ::    r4_value =   (/2.71828/)
    real*8, parameter, dimension(*) ::    r8_value =   (/3.141592653589793238462643383279502884197169/)

    integer*8 :: file_id, status

#if defined(PARALLEL_IO)
    integer :: ierr
    call mpi_init(ierr)
#endif

    ! abort program on any H5hut error
    call h5_abort_on_error()

    call h5_set_verbosity_level (verbosity_level)

    !  MPI_COMM_WORLD is used, if file is opened with default properties
    file_id = h5_openfile (FNAME, H5_O_WRONLY, H5_PROP_DEFAULT)

    ! open step 1
    status = h5_setstep (file_id, int8(1))

    ! write attributes
    status = h5_writestepattrib_string (file_id, ATTR_STRING, string_value)
    status = h5_writestepattrib_i4     (file_id, ATTR_I4,     i4_value, int8(size(i4_value, 1)))
    status = h5_writestepattrib_i8     (file_id, ATTR_I8,     i8_value, int8(size(i8_value, 1)))
    status = h5_writestepattrib_r4     (file_id, ATTR_R4,     r4_value, int8(size(r4_value, 1)))
    status = h5_writestepattrib_r8     (file_id, ATTR_R8,     r8_value, int8(size(r8_value, 1)))

    ! cleanup
    status = h5_closefile (file_id)

#if defined(PARALLEL_IO)
    call mpi_finalize(ierr)
#endif
  end program write_step_attribs


text:
c*************************************************************************
c                         TU4_VKICKB.F
c*************************************************************************
c VKICKB gives a 'kick' to the bary. velocities for use in the symplectic
c code based on Ham. of the form T(p) + V(x).
c
c             Input:
c                 nbod          ==>  number of massive bodies (int scalar)
c                 vxb,vyb,vzb   ==>  initial velocity in beri coord
c                                    (real arrays)
c                 axb,ayb,azb   ==>  accel in beri coord (real arrays)
c                 dt            ==>  time step
c             Output:
c                 vxb,vyb,vzb   ==>  final velocity in beri coord
c                                       (real arrays)
c
c Remarks:  Based on Martin's NB4M routines
c Authors:  Martin Duncan
c Date:    3/8/93
c Last revision:

	  subroutine tu4_vkickb(nbod,vxb,vyb,vzb,axb,ayb,azb,dt)

      include '../swift.inc'

c...  Inputs Only:
      integer nbod
      real*8 axb(nbod),ayb(nbod),azb(nbod),dt

c...  Inputs and Outputs:
      real*8 vxb(nbod),vyb(nbod),vzb(nbod)

c...  Internals
      integer i

c---
c...  Executable code

	  do i=1,nbod
	     vxb(i) = vxb(i) + axb(i)*dt
	     vyb(i) = vyb(i) + ayb(i)*dt
	     vzb(i) = vzb(i) + azb(i)*dt
	  enddo

	  return
	  end   ! tu4_vkickb
c____________________________________________________________________________

max_stars_repo_name: Simske/exostriker