text:
C
C**************************  SHLPIC ************************************
C
C  restore diagonal elements of LHS system matrix corresponding to
C  Dirichlet nodes.
C  Also, set solution at Dirichlet nodes to the prescribed
C  values. This is done since the solved solution at
C  Dirichlet nodes may not be exactly equal to the presribed
C  values, due to inaccuracies and roundoff errors which could
C  arise from the way we imposed Dirichlet conditions (by multiplying
C  diagonal terms by a 'large' number).
C  Picard scheme (symmetric storage).
C
C***********************************************************************
C
      SUBROUTINE SHLPIC(NP,TOPOL,CONTP,PRESC,LHSP,COEF1,PNEW,POLD,
     1                  NSF,NSFNUM,NSFNOD,SFEX,LHSSF,DUPUIT,Z,
     2                  NNOD,IFATM,LHSATM)
C
      IMPLICIT  NONE
      INCLUDE  'CATHY.H'
      INTEGER   I,J,IND,INOD
      INTEGER   NP,NSF,NNOD,DUPUIT
      INTEGER   TOPOL(*),CONTP(*),NSFNUM(*),NSFNOD(NSFMAX,*)
      INTEGER   SFEX(NSFMAX,*)
      INTEGER   IFATM(*)
      REAL*8    PRESC(*),LHSP(*),COEF1(*),PNEW(*),POLD(*),Z(*)
      REAL*8    LHSSF(NSFMAX,*),LHSATM(*)
C
      DO I=1,NP
         J=CONTP(I)
         PNEW(J)=PRESC(I)
         IND=TOPOL(J)
         COEF1(IND)=LHSP(I)
      END DO
      DO I=1,NSF
         DO J=1,NSFNUM(I)
            IF (SFEX(I,J).EQ.1) THEN
               INOD=NSFNOD(I,J)
               IND=TOPOL(INOD)
c               IF (DUPUIT.EQ.0) THEN
                  PNEW(INOD)=0.0D0
c               ELSE
c            PNEW(INOD)=0.0D0+Z(NSFNOD(I,SFEX(I)))-Z(INOD)
c               END IF
               COEF1(IND)=LHSSF(I,J)
            END IF
         END DO
      END DO
      DO I=1,NNOD
         IF (IFATM(I) .EQ. 1 .OR. IFATM(I) .EQ. 2) THEN
            IND=TOPOL(I)
            PNEW(I)=POLD(I)
            COEF1(IND)=LHSATM(I)
         END IF
      END DO
C
C     WRITE(66666,*)
      RETURN
      END


!simple fortran 95 program written from GNU nano, lets take a trip back to the 60's aye :D

program main
    implicit none
    CHARACTER(LEN=15) :: name
    integer :: age
    !could not induce you're as it would cause a syntax error
    print *, '==== Enter your name ===='
    read(*,*) name
    print *, '==== Enter your age  ===='
    print *, " > "
    read(*,*) age
    print *, 'Hello there ', name
    IF (age > 18) THEN
       print *, 'age => ', age
       print *, 'you are NOT below 18'
    ELSE IF (age < 18) THEN
       print *, 'age => ', age
       print *, 'you are blow 18'
 END IF

end program main


text:
      SUBROUTINE HOSETG( NOUT, LERR, SHORTO,
     &             SHYDRA, FHYDRA, LLINE,
     &             SPDBSP, FPDBSP, LSPHBX,
     &             MOLQPT, ATMAX, ATNO, ATXYZ, ATBND)
      IMPLICIT NONE
C ********************************************************************
C *                                                                  *
C * This software is an unpublished work containing confidential and *
C * proprietary information of Birkbeck College. Use, disclosure,    *
C * reproduction and transfer of this work without the express       *
C * written consent of Birkbeck College are prohibited. This notice  *
C * must be attached to all copies or extracts of the software.      *
C *                                                                  *
C * (c) 1995 Oliver Smart & Birkbeck College, All rights reserved    *
C * (c) 1996 Oliver Smart & Birkbeck College, All rights reserved    *
C *                                                                  *
C ********************************************************************
C
C Modification history:
C
C Date  Author          Modification
C 11/95 O.S. Smart      First version
C 01/00 O.S.		S/r NEWOP introduced
C
C This s/r deals with opening the graphics files at the start of
C a HOLE run.

C output stream no (to user)
C (return unchanged)
      INTEGER                   NOUT

C if an error found is found in s/r - set lerr true and program will stop,
      LOGICAL                   LERR

C  variable to reduce text output -
C if oshort =
C 0: full text output
C 1: no in run progress reported
C 2: ditto plus no graph type output - only leaving
C    Minimum radius & conductance calculations.
C 3: All text output other than input card mirroring
C    and error messages turned off
C (return unchanged)
      INTEGER                   SHORTO

C quanta 3D plot format (.qpt) file stream, name
      INTEGER			SHYDRA
      CHARACTER*200		FHYDRA
C logical flag to indicate whether we are drawing lines from
C the HOLE sphere surfaces to closest atom
C (return unchanged)
      LOGICAL			LLINE

C sphere centre pdb format file (for use with sphqpb)
C stream number, filename
      INTEGER			SPDBSP
      CHARACTER*200		FPDBSP

C sphere-box option on?
C (return unchanged)
      LOGICAL			LSPHBX

C Output filename for .qpt format stick plot of the molecule
C (return unchanged)
      CHARACTER*200             MOLQPT

C maximum number of atoms
C (return unchanged)
      INTEGER                   ATMAX

C the actual number of atoms
C (return unchanged)
      INTEGER                   ATNO

C atom's Cartesian coords
C (return unchanged)
      DOUBLE PRECISION          ATXYZ(3,ATMAX)

C atom's bonding radii - if two atoms are closer than
C the sum of their bonding radii they are regarded as bonded
C and a line is drawn
C (return unchanged)
      DOUBLE PRECISION          ATBND(ATMAX)

C internal vbles ***************

C end of decs ******************

C open the plot file - if we need to
      IF (FHYDRA(1:4).NE.'NONE') THEN
C open routine allows the overwriting of existing files
C (last arguement set true for binary files)
        CALL NEWOP( NOUT, LERR, SHORTO,
     &                  FHYDRA, SHYDRA, .TRUE.)
        IF (LERR) THEN
C error cannot open hydra binary plot output file
        WRITE( NOUT,  '(A)')
     &' ERROR',
     &' Cannot open hydra binary plot output file:  '//
     &   FHYDRA(1:INDEX( FHYDRA,'      ')-1),
     &' Aborting'
         GOTO 55555
        ENDIF

C tell user:
        INQUIRE( SHYDRA, NAME = FHYDRA)
        IF (SHORTO.LT.3) WRITE( NOUT,  '(/1X,A)')
     &    'Have opened hydra binary plot output file:  '//
     &    FHYDRA(1:INDEX( FHYDRA,'      ')-1)

C colour 1 is for lines joining sphere edge to closest atoms
C Change 26/10/93 do in colour 10 with alternative qplot 20
        IF (LLINE) WRITE(SHYDRA) 1.0, 10.0, -55.0, 20.0

      ELSE
        WRITE(NOUT,*) 'Will not produce graphical output'
      ENDIF

C open the pdb format file for sphere centres - if we need to
      IF (FPDBSP(1:4).NE.'NONE') THEN
C open routine allows the overwriting of existing files
C (last arguement set true for binary files)
        CALL NEWOP( NOUT, LERR, SHORTO,
     &                  FPDBSP, SPDBSP, .FALSE.)
        IF (LERR) THEN
C error cannot open pdb sphere centre output file
          WRITE( NOUT,  '(A)')
     &  ' ERROR',
     &  ' cannot pdb sphere centre file:  '//
     &   FPDBSP(1:INDEX( FPDBSP,'      ')-1)
          GOTO 55555
        ENDIF
C tell user output name
        INQUIRE( SPDBSP, NAME = FPDBSP)
        IF (SHORTO.LT.3) WRITE( NOUT,  '(/1X,A)')
     &'Have opened file for pdb format output of sphere centres:  '//
     &    FPDBSP(1:INDEX( FPDBSP,'      ')-1)
C if doing sphere-box record the fact
        IF (LSPHBX) WRITE( SPDBSP, '(A)') 'SPHERE-BOX-ON'

      ELSE
        IF (SHORTO.LT.3) WRITE( NOUT, '(A)')
     &' Will not produce file with sphere centre',
     &'      information in pdb format'
      ENDIF

C 30 October 1995 molqpt can be specified to write
C a .qpt file of original molecule
      IF (MOLQPT(1:4).NE.'NONE') THEN
        CALL WMOLQP( NOUT, LERR, MOLQPT, ATMAX, ATNO, ATXYZ, ATBND,
     &		     SHORTO)
        IF (LERR) GOTO 55555
      ENDIF

55555 RETURN

      END

ext: f


text:
great_rsd(u1,kk1).
great_rsd(aa1,c1).
great_rsd(w1,x1).
great_rsd(f1,v1).
great_rsd(a1,v1).
great_rsd(a1,t1).
great_rsd(w1,ff1).
great_rsd(a1,m1).
great_rsd(u1,k1).
great_rsd(p1,x1).
great_rsd(aa1,q1).
great_rsd(z1,jj1).
great_rsd(k1,jj1).
great_rsd(k1,dd1).
great_rsd(a1,x1).
great_rsd(u1,hh1).
great_rsd(q1,ii1).
great_rsd(aa1,m1).
great_rsd(e1,jj1).
great_rsd(h1,x1).
great_rsd(ff1,g1).
great_rsd(aa1,dd1).
great_rsd(g1,h1).
great_rsd(aa1,k1).
great_rsd(aa1,l1).
great_rsd(ee1,jj1).
great_rsd(z1,h1).
great_rsd(kk1,z1).
great_rsd(k1,hh1).
great_rsd(f1,a1).
great_rsd(a1,u1).
great_rsd(a1,jj1).
great_rsd(cc1,h1).
great_rsd(aa1,h1).
great_rsd(ee1,h1).
great_rsd(g1,x1).
great_rsd(w1,t1).
great_rsd(m1,dd1).
great_rsd(e1,x1).
great_rsd(w1,c1).
great_rsd(kk1,x1).
great_rsd(m1,g1).
great_rsd(a1,c1).
great_rsd(m1,hh1).
great_rsd(a1,k1).
great_rsd(c1,dd1).
great_rsd(g1,l1).
great_rsd(e1,q1).
great_rsd(ff1,z1).
great_rsd(ff1,c1).
great_rsd(w1,hh1).
great_rsd(a1,h1).
great_rsd(e1,c1).
great_rsd(v1,c1).
great_rsd(g1,ii1).
great_rsd(f1,g1).
great_rsd(v1,x1).
great_rsd(l1,x1).
great_rsd(p1,l1).
great_rsd(v1,m1).
great_rsd(m1,z1).
great_rsd(w1,m1).
great_rsd(e1,v1).
great_rsd(c1,ii1).
great_rsd(q1,ff1).
great_rsd(v1,q1).
great_rsd(cc1,ff1).
great_rsd(m1,kk1).
great_rsd(v1,g1).
great_rsd(e1,p1).
great_rsd(aa1,g1).
great_rsd(ff1,jj1).
great_rsd(q1,hh1).
great_rsd(aa1,ee1).
great_rsd(f1,ii1).
great_rsd(q1,m1).
great_rsd(ff1,l1).
great_rsd(t1,m1).
great_rsd(u1,ee1).
great_rsd(k1,c1).
great_rsd(h1,jj1).
great_rsd(f1,h1).
great_rsd(e1,g1).
great_rsd(a1,ff1).
great_rsd(p1,u1).
great_rsd(p1,w1).
great_rsd(z1,x1).
great_rsd(cc1,dd1).
great_rsd(e1,w1).
great_rsd(ii1,h1).
great_rsd(u1,h1).
great_rsd(u1,dd1).
great_rsd(cc1,u1).
great_rsd(v1,z1).
great_rsd(v1,kk1).
great_rsd(u1,z1).
great_rsd(m1,c1).
great_rsd(ii1,x1).
great_rsd(p1,m1).
great_rsd(t1,x1).
great_rsd(p1,aa1).
great_rsd(cc1,x1).
great_rsd(a1,hh1).
great_rsd(w1,q1).
great_rsd(e1,k1).
great_rsd(a1,cc1).
great_rsd(f1,l1).
great_rsd(a1,w1).
great_rsd(q1,g1).
great_rsd(t1,aa1).
great_rsd(f1,u1).
great_rsd(p1,g1).
great_rsd(m1,l1).
great_rsd(a1,dd1).
great_rsd(t1,hh1).
great_rsd(f1,x1).
great_rsd(k1,l1).
great_rsd(cc1,kk1).
great_rsd(w1,h1).
great_rsd(ii1,l1).
great_rsd(ff1,m1).
great_rsd(g1,c1).
great_rsd(ee1,x1).
great_rsd(v1,ff1).
great_rsd(p1,a1).
great_rsd(v1,aa1).
great_rsd(v1,k1).
great_rsd(e1,ff1).
great_rsd(p1,kk1).
great_rsd(p1,jj1).
great_rsd(q1,jj1).
great_rsd(c1,l1).
great_rsd(cc1,hh1).
great_rsd(a1,q1).
great_rsd(f1,p1).
great_rsd(l1,dd1).
great_rsd(u1,x1).
great_rsd(m1,ii1).
great_rsd(cc1,z1).
great_rsd(p1,q1).
great_rsd(t1,k1).
great_rsd(f1,k1).
great_rsd(u1,g1).
great_rsd(z1,dd1).
great_rsd(t1,q1).
great_rsd(k1,x1).
great_rsd(k1,ii1).
great_rsd(q1,h1).
great_rsd(f1,ee1).
great_rsd(aa1,z1).
great_rsd(c1,jj1).
great_rsd(p1,c1).
great_rsd(cc1,q1).
great_rsd(t1,z1).
great_rsd(e1,z1).
great_rsd(cc1,t1).
great_rsd(f1,w1).
great_rsd(t1,c1).
great_rsd(f1,dd1).
great_rsd(p1,dd1).
great_rsd(kk1,c1).
great_rsd(p1,v1).
great_rsd(w1,l1).
great_rsd(w1,ee1).
great_rsd(t1,l1).
great_rsd(m1,jj1).
great_rsd(u1,jj1).
great_rsd(cc1,c1).
great_rsd(t1,jj1).
great_rsd(a1,l1).
great_rsd(w1,jj1).
great_rsd(hh1,jj1).
great_rsd(a1,kk1).
great_rsd(q1,dd1).
great_rsd(q1,z1).
great_rsd(t1,g1).
great_rsd(ee1,l1).
great_rsd(cc1,v1).
great_rsd(a1,ee1).
great_rsd(jj1,x1).
great_rsd(ee1,hh1).
great_rsd(c1,h1).
great_rsd(f1,m1).
great_rsd(u1,aa1).
great_rsd(ii1,dd1).
great_rsd(u1,t1).
great_rsd(cc1,l1).
great_rsd(ff1,k1).
great_rsd(v1,h1).
great_rsd(u1,ff1).
great_rsd(h1,l1).
great_rsd(aa1,hh1).
great_rsd(ff1,kk1).
great_rsd(a1,ii1).
great_rsd(v1,ii1).
great_rsd(t1,ff1).
great_rsd(h1,dd1).
great_rsd(q1,x1).
great_rsd(ff1,x1).
great_rsd(e1,t1).
great_rsd(f1,hh1).
great_rsd(kk1,l1).
great_rsd(ee1,dd1).
great_rsd(e1,aa1).
great_rsd(v1,u1).
great_rsd(f1,ff1).
great_rsd(c1,z1).
great_rsd(aa1,x1).
great_rsd(p1,k1).
great_rsd(g1,hh1).
great_rsd(dd1,x1).
great_rsd(m1,x1).
great_rsd(p1,hh1).
great_rsd(g1,z1).
great_rsd(w1,k1).
great_rsd(ii1,k1).
great_rsd(m1,aa1).
great_rsd(w1,cc1).
great_rsd(jj1,kk1).
great_rsd(x1,ii1).
great_rsd(g1,t1).
great_rsd(m1,w1).
great_rsd(u1,w1).
great_rsd(aa1,t1).
great_rsd(c1,p1).
great_rsd(v1,e1).
great_rsd(k1,ee1).
great_rsd(z1,ee1).
great_rsd(l1,p1).
great_rsd(ff1,q1).
great_rsd(ee1,a1).
great_rsd(g1,m1).
great_rsd(x1,m1).
great_rsd(x1,v1).
great_rsd(h1,p1).
great_rsd(aa1,p1).
great_rsd(dd1,cc1).
great_rsd(cc1,e1).
great_rsd(hh1,f1).
great_rsd(z1,c1).
great_rsd(ff1,t1).
great_rsd(ii1,ff1).
great_rsd(q1,a1).
great_rsd(dd1,ee1).
great_rsd(c1,aa1).
great_rsd(ii1,a1).
great_rsd(l1,h1).
great_rsd(v1,a1).
great_rsd(jj1,u1).
great_rsd(hh1,a1).
great_rsd(h1,f1).
great_rsd(dd1,l1).
great_rsd(ee1,cc1).
great_rsd(ff1,ee1).
great_rsd(w1,p1).
great_rsd(dd1,h1).
great_rsd(jj1,m1).
great_rsd(c1,kk1).
great_rsd(k1,u1).
great_rsd(kk1,aa1).
great_rsd(dd1,jj1).
great_rsd(ee1,t1).
great_rsd(q1,cc1).
great_rsd(g1,f1).
great_rsd(ii1,ee1).
great_rsd(z1,f1).
great_rsd(jj1,ee1).
great_rsd(g1,ee1).
great_rsd(h1,c1).
great_rsd(x1,h1).
great_rsd(ii1,p1).
great_rsd(v1,p1).
great_rsd(u1,cc1).
great_rsd(dd1,aa1).
great_rsd(k1,f1).
great_rsd(q1,v1).
great_rsd(z1,m1).
great_rsd(jj1,k1).
great_rsd(q1,e1).
great_rsd(m1,a1).
great_rsd(jj1,v1).
great_rsd(u1,p1).
great_rsd(h1,hh1).
great_rsd(hh1,q1).
great_rsd(g1,q1).
great_rsd(ii1,kk1).
great_rsd(h1,v1).
great_rsd(dd1,w1).
great_rsd(k1,a1).
great_rsd(ff1,u1).
great_rsd(z1,p1).
great_rsd(k1,q1).
great_rsd(c1,g1).
great_rsd(v1,f1).
great_rsd(jj1,ff1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
This is a good one for those lazy Tuesday evenings.. You know, the kind where you dont have anything to do, and Gone in 60 seconds is playing on TV for the 40billionth time. Con Air, Face/Off or The Rock will also work.

First, pour yourself as many shots as you have glasses for. This could get very ugly very quickly if you use hard liquor, so Id suggest starting with beer or some sort of mixed shot.

Second, start watching the movie.

Now, heres the important part: Watch Nicholas Cage closely. Every time he looks perplexed, TAKE A SHOT.

Just in case youre a visual learner, heres a diagram to help you:

Whoever makes it the longest wins! (its pretty tough to make it through the whole film)

Or, you can play alone. When you pass out, you win! But its less pathetic with at least one friend.

.

.

.

(btw, read the warning at the top of Drinking games. Nicholas Cage doesnt like it when you get alcohol poisoning..)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine echarge1  --  charge-charge energy & derivs  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "echarge1" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c
c
      subroutine echarge1
      implicit none
      include 'sizes.i'
      include 'cutoff.i'
      include 'warp.i'
c
c
c     choose the method for summing over pairwise interactions
c
      if (use_smooth) then
         call echarge1g
      else if (use_ewald) then
         if (use_clist) then
            call echarge1f
         else if (use_lights) then
            call echarge1e
         else
            call echarge1d
         end if
      else if (use_clist) then
         call echarge1c
      else if (use_lights) then
         call echarge1b
      else
         call echarge1a
      end if
      return
      end
c
c
c     ###############################################################
c     ##                                                           ##
c     ##  subroutine echarge1a  --  charge derivs via double loop  ##
c     ##                                                           ##
c     ###############################################################
c
c
c     "echarge1a" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     using a pairwise double loop
c
c
      subroutine echarge1a
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'cell.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'group.i'
      include 'inter.i'
      include 'molcul.i'
      include 'shunt.i'
      include 'usage.i'
      include 'virial.i'
      integer i,j,k
      integer ii,kk
      integer in,kn
      integer ic,kc
      real*8 e,fgrp,de,dc
      real*8 f,fi,fik
      real*8 r,r2,rb,rb2
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 xc,yc,zc
      real*8 xic,yic,zic
      real*8 dedx,dedy,dedz
      real*8 dedxc,dedyc,dedzc
      real*8 shift,taper,trans
      real*8 dtaper,dtrans
      real*8 rc,rc2,rc3,rc4
      real*8 rc5,rc6,rc7
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      real*8, allocatable :: cscale(:)
      logical proceed,usei
      character*6 mode
c
c
c     zero out the charge interaction energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     set conversion factor, cutoff and switching coefficients
c
      f = electric / dielec
      mode = 'CHARGE'
      call switch (mode)
c
c     compute the charge interaction energy and first derivatives
c
      do ii = 1, nion-1
         i = iion(ii)
         in = jion(ii)
         ic = kion(ii)
         xic = x(ic)
         yic = y(ic)
         zic = z(ic)
         xi = x(i) - xic
         yi = y(i) - yic
         zi = z(i) - zic
         fi = f * pchg(ii)
         usei = (use(i) .or. use(ic))
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kk = ii+1, nion
            k = iion(kk)
            kn = jion(kk)
            kc = kion(kk)
            proceed = .true.
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            if (proceed)  proceed = (usei .or. use(k) .or. use(kc))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xc = xic - x(kc)
               yc = yic - y(kc)
               zc = zic - z(kc)
               if (use_bounds)  call image (xc,yc,zc)
               rc2 = xc*xc + yc*yc + zc*zc
               if (rc2 .le. off2) then
                  xr = xc + xi - x(k) + x(kc)
                  yr = yc + yi - y(k) + y(kc)
                  zr = zc + zi - z(k) + z(kc)
                  r2 = xr*xr + yr*yr + zr*zr
                  r = sqrt(r2)
                  rb = r + ebuffer
                  rb2 = rb * rb
                  fik = fi * pchg(kk) * cscale(kn)
                  e = fik / rb
                  de = -fik / rb2
                  dc = 0.0d0
c
c     use shifted energy switching if near the cutoff distance
c
                  shift = fik / (0.5d0*(off+cut))
                  e = e - shift
                  if (rc2 .gt. cut2) then
                     rc = sqrt(rc2)
                     rc3 = rc2 * rc
                     rc4 = rc2 * rc2
                     rc5 = rc2 * rc3
                     rc6 = rc3 * rc3
                     rc7 = rc3 * rc4
                     taper = c5*rc5 + c4*rc4 + c3*rc3
     &                          + c2*rc2 + c1*rc + c0
                     dtaper = 5.0d0*c5*rc4 + 4.0d0*c4*rc3
     &                           + 3.0d0*c3*rc2 + 2.0d0*c2*rc + c1
                     trans = fik * (f7*rc7 + f6*rc6 + f5*rc5 + f4*rc4
     &                               + f3*rc3 + f2*rc2 + f1*rc + f0)
                     dtrans = fik * (7.0d0*f7*rc6 + 6.0d0*f6*rc5
     &                               + 5.0d0*f5*rc4 + 4.0d0*f4*rc3
     &                             + 3.0d0*f3*rc2 + 2.0d0*f2*rc + f1)
                     dc = (e*dtaper + dtrans) / rc
                     de = de * taper
                     e = e*taper + trans
                  end if
c
c     scale the interaction based on its group membership
c
                  if (use_group) then
                     e = e * fgrp
                     de = de * fgrp
                     dc = dc * fgrp
                  end if
c
c     form the chain rule terms for derivative expressions
c
                  de = de / r
                  dedx = de * xr
                  dedy = de * yr
                  dedz = de * zr
                  dedxc = dc * xc
                  dedyc = dc * yc
                  dedzc = dc * zc
c
c     increment the overall energy and derivative expressions
c
                  ec = ec + e
                  dec(1,i) = dec(1,i) + dedx
                  dec(2,i) = dec(2,i) + dedy
                  dec(3,i) = dec(3,i) + dedz
                  dec(1,ic) = dec(1,ic) + dedxc
                  dec(2,ic) = dec(2,ic) + dedyc
                  dec(3,ic) = dec(3,ic) + dedzc
                  dec(1,k) = dec(1,k) - dedx
                  dec(2,k) = dec(2,k) - dedy
                  dec(3,k) = dec(3,k) - dedz
                  dec(1,kc) = dec(1,kc) - dedxc
                  dec(2,kc) = dec(2,kc) - dedyc
                  dec(3,kc) = dec(3,kc) - dedzc
c
c     increment the internal virial tensor components
c
                  vxx = xr*dedx + xc*dedxc
                  vyx = yr*dedx + yc*dedxc
                  vzx = zr*dedx + zc*dedxc
                  vyy = yr*dedy + yc*dedyc
                  vzy = zr*dedy + zc*dedyc
                  vzz = zr*dedz + zc*dedzc
                  vir(1,1) = vir(1,1) + vxx
                  vir(2,1) = vir(2,1) + vyx
                  vir(3,1) = vir(3,1) + vzx
                  vir(1,2) = vir(1,2) + vyx
                  vir(2,2) = vir(2,2) + vyy
                  vir(3,2) = vir(3,2) + vzy
                  vir(1,3) = vir(1,3) + vzx
                  vir(2,3) = vir(2,3) + vzy
                  vir(3,3) = vir(3,3) + vzz
c
c     increment the total intermolecular energy
c
                  if (molcule(i) .ne. molcule(k)) then
                     einter = einter + e
                  end if
               end if
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     for periodic boundary conditions with large cutoffs
c     neighbors must be found by the replicates method
c
      if (.not. use_replica)  return
c
c     calculate interaction energy with other unit cells
c
      do ii = 1, nion
         i = iion(ii)
         in = jion(ii)
         ic = kion(ii)
         usei = (use(i) .or. use(ic))
         xic = x(ic)
         yic = y(ic)
         zic = z(ic)
         xi = x(i) - xic
         yi = y(i) - yic
         zi = z(i) - zic
         fi = f * pchg(ii)
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kk = ii, nion
            k = iion(kk)
            kn = jion(kk)
            kc = kion(kk)
            proceed = .true.
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            if (proceed)  proceed = (usei .or. use(k) .or. use(kc))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               do j = 1, ncell
                  xc = xic - x(kc)
                  yc = yic - y(kc)
                  zc = zic - z(kc)
                  call imager (xc,yc,zc,j)
                  rc2 = xc*xc + yc*yc + zc*zc
                  if (rc2 .le. off2) then
                     xr = xc + xi - x(k) + x(kc)
                     yr = yc + yi - y(k) + y(kc)
                     zr = zc + zi - z(k) + z(kc)
                     r2 = xr*xr + yr*yr + zr*zr
                     r = sqrt(r2)
                     rb = r + ebuffer
                     rb2 = rb * rb
                     fik = fi * pchg(kk)
                     if (use_polymer) then
                        if (r2 .le. polycut2)  fik = fik * cscale(kn)
                     end if
                     e = fik / rb
                     de = -fik / rb2
                     dc = 0.0d0
c
c     use shifted energy switching if near the cutoff distance
c
                     shift = fik / (0.5d0*(off+cut))
                     e = e - shift
                     if (rc2 .gt. cut2) then
                        rc = sqrt(rc2)
                        rc3 = rc2 * rc
                        rc4 = rc2 * rc2
                        rc5 = rc2 * rc3
                        rc6 = rc3 * rc3
                        rc7 = rc3 * rc4
                        taper = c5*rc5 + c4*rc4 + c3*rc3
     &                             + c2*rc2 + c1*rc + c0
                        dtaper = 5.0d0*c5*rc4 + 4.0d0*c4*rc3
     &                              + 3.0d0*c3*rc2 + 2.0d0*c2*rc + c1
                        trans = fik * (f7*rc7 + f6*rc6 + f5*rc5 + f4*rc4
     &                                  + f3*rc3 + f2*rc2 + f1*rc + f0)
                        dtrans = fik * (7.0d0*f7*rc6 + 6.0d0*f6*rc5
     &                                  + 5.0d0*f5*rc4 + 4.0d0*f4*rc3
     &                                + 3.0d0*f3*rc2 + 2.0d0*f2*rc + f1)
                        dc = (e*dtaper + dtrans) / rc
                        de = de * taper
                        e = e*taper + trans
                     end if
c
c     scale the interaction based on its group membership
c
                     if (use_group) then
                        e = e * fgrp
                        de = de * fgrp
                        dc = dc * fgrp
                     end if
c
c     form the chain rule terms for derivative expressions
c
                     de = de / r
                     dedx = de * xr
                     dedy = de * yr
                     dedz = de * zr
                     dedxc = dc * xc
                     dedyc = dc * yc
                     dedzc = dc * zc
c
c     increment the energy and gradient values
c
                     if (i .eq. k)  e = 0.5d0 * e
                     ec = ec + e
                     dec(1,i) = dec(1,i) + dedx
                     dec(2,i) = dec(2,i) + dedy
                     dec(3,i) = dec(3,i) + dedz
                     dec(1,ic) = dec(1,ic) + dedxc
                     dec(2,ic) = dec(2,ic) + dedyc
                     dec(3,ic) = dec(3,ic) + dedzc
                     if (i .ne. k) then
                        dec(1,k) = dec(1,k) - dedx
                        dec(2,k) = dec(2,k) - dedy
                        dec(3,k) = dec(3,k) - dedz
                        dec(1,kc) = dec(1,kc) - dedxc
                        dec(2,kc) = dec(2,kc) - dedyc
                        dec(3,kc) = dec(3,kc) - dedzc
                     end if
c
c     increment the internal virial tensor components
c
                     vxx = xr*dedx + xc*dedxc
                     vyx = yr*dedx + yc*dedxc
                     vzx = zr*dedx + zc*dedxc
                     vyy = yr*dedy + yc*dedyc
                     vzy = zr*dedy + zc*dedyc
                     vzz = zr*dedz + zc*dedzc
                     vir(1,1) = vir(1,1) + vxx
                     vir(2,1) = vir(2,1) + vyx
                     vir(3,1) = vir(3,1) + vzx
                     vir(1,2) = vir(1,2) + vyx
                     vir(2,2) = vir(2,2) + vyy
                     vir(3,2) = vir(3,2) + vzy
                     vir(1,3) = vir(1,3) + vzx
                     vir(2,3) = vir(2,3) + vzy
                     vir(3,3) = vir(3,3) + vzz
c
c     increment the total intermolecular energy
c
                     einter = einter + e
                  end if
               end do
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine echarge1b  --  method of lights charge derivs  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "echarge1b" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     using the method of lights
c
c
      subroutine echarge1b
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'boxes.i'
      include 'cell.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'group.i'
      include 'inter.i'
      include 'light.i'
      include 'molcul.i'
      include 'shunt.i'
      include 'usage.i'
      include 'virial.i'
      integer i,j,k,ii,kk
      integer in,ic,kn,kc
      integer kgy,kgz,kmap
      integer start,stop
      real*8 e,fgrp,de,dc
      real*8 f,fi,fik
      real*8 r,r2,rb,rb2
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 xc,yc,zc
      real*8 xic,yic,zic
      real*8 dedx,dedy,dedz
      real*8 dedxc,dedyc,dedzc
      real*8 shift,taper,trans
      real*8 dtaper,dtrans
      real*8 rc,rc2,rc3,rc4
      real*8 rc5,rc6,rc7
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      real*8, allocatable :: cscale(:)
      real*8, allocatable :: xsort(:)
      real*8, allocatable :: ysort(:)
      real*8, allocatable :: zsort(:)
      logical proceed,usei
      logical prime,repeat
      character*6 mode
c
c
c     zero out the charge interaction energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
      allocate (xsort(8*n))
      allocate (ysort(8*n))
      allocate (zsort(8*n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     set conversion factor, cutoff and switching coefficients
c
      f = electric / dielec
      mode = 'CHARGE'
      call switch (mode)
c
c     transfer the interaction site coordinates to sorting arrays
c
      do i = 1, nion
         k = kion(i)
         xsort(i) = x(k)
         ysort(i) = y(k)
         zsort(i) = z(k)
      end do
c
c     use the method of lights to generate neighbors
c
      call lights (off,nion,xsort,ysort,zsort)
c
c     loop over all atoms computing the interactions
c
      do ii = 1, nion
         i = iion(ii)
         in = jion(ii)
         ic = kion(ii)
         xic = xsort(rgx(ii))
         yic = ysort(rgy(ii))
         zic = zsort(rgz(ii))
         xi = x(i) - x(ic)
         yi = y(i) - y(ic)
         zi = z(i) - z(ic)
         fi = f * pchg(ii)
         usei = (use(i) .or. use(ic))
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     loop over method of lights neighbors of current atom
c
         if (kbx(ii) .le. kex(ii)) then
            repeat = .false.
            start = kbx(ii) + 1
            stop = kex(ii)
         else
            repeat = .true.
            start = 1
            stop = kex(ii)
         end if
   10    continue
         do j = start, stop
            kk = locx(j)
            kgy = rgy(kk)
            if (kby(ii) .le. key(ii)) then
               if (kgy.lt.kby(ii) .or. kgy.gt.key(ii))  goto 20
            else
               if (kgy.lt.kby(ii) .and. kgy.gt.key(ii))  goto 20
            end if
            kgz = rgz(kk)
            if (kbz(ii) .le. kez(ii)) then
               if (kgz.lt.kbz(ii) .or. kgz.gt.kez(ii))  goto 20
            else
               if (kgz.lt.kbz(ii) .and. kgz.gt.kez(ii))  goto 20
            end if
            kmap = kk - ((kk-1)/nion)*nion
            k = iion(kmap)
            kn = jion(kmap)
            kc = kion(kmap)
            prime = (kk .le. nion)
c
c     decide whether to compute the current interaction
c
            proceed = .true.
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            if (proceed)  proceed = (usei .or. use(k) .or. use(kc))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xc = xic - xsort(j)
               yc = yic - ysort(kgy)
               zc = zic - zsort(kgz)
               if (use_bounds) then
                  if (abs(xc) .gt. xcell2)  xc = xc - sign(xcell,xc)
                  if (abs(yc) .gt. ycell2)  yc = yc - sign(ycell,yc)
                  if (abs(zc) .gt. zcell2)  zc = zc - sign(zcell,zc)
                  if (monoclinic) then
                     xc = xc + zc*beta_cos
                     zc = zc * beta_sin
                  else if (triclinic) then
                     xc = xc + yc*gamma_cos + zc*beta_cos
                     yc = yc*gamma_sin + zc*beta_term
                     zc = zc * gamma_term
                  end if
               end if
               rc2 = xc*xc + yc*yc + zc*zc
               if (rc2 .le. off2) then
                  xr = xc + xi - x(k) + x(kc)
                  yr = yc + yi - y(k) + y(kc)
                  zr = zc + zi - z(k) + z(kc)
                  r2 = xr*xr + yr*yr + zr*zr
                  r = sqrt(r2)
                  rb = r + ebuffer
                  rb2 = rb * rb
                  fik = fi * pchg(kmap)
                  if (prime)  fik = fik * cscale(kn)
                  if (use_polymer) then
                     if (r2 .gt. polycut2)  fik = fi * pchg(kmap)
                  end if
                  e = fik / rb
                  de = -fik / rb2
                  dc = 0.0d0
c
c     use shifted energy switching if near the cutoff distance
c
                  shift = fik / (0.5d0*(off+cut))
                  e = e - shift
                  if (rc2 .gt. cut2) then
                     rc = sqrt(rc2)
                     rc3 = rc2 * rc
                     rc4 = rc2 * rc2
                     rc5 = rc2 * rc3
                     rc6 = rc3 * rc3
                     rc7 = rc3 * rc4
                     taper = c5*rc5 + c4*rc4 + c3*rc3
     &                          + c2*rc2 + c1*rc + c0
                     dtaper = 5.0d0*c5*rc4 + 4.0d0*c4*rc3
     &                           + 3.0d0*c3*rc2 + 2.0d0*c2*rc + c1
                     trans = fik * (f7*rc7 + f6*rc6 + f5*rc5 + f4*rc4
     &                               + f3*rc3 + f2*rc2 + f1*rc + f0)
                     dtrans = fik * (7.0d0*f7*rc6 + 6.0d0*f6*rc5
     &                               + 5.0d0*f5*rc4 + 4.0d0*f4*rc3
     &                             + 3.0d0*f3*rc2 + 2.0d0*f2*rc + f1)
                     dc = (e*dtaper + dtrans) / rc
                     de = de * taper
                     e = e*taper + trans
                  end if
c
c     scale the interaction based on its group membership
c
                  if (use_group) then
                     e = e * fgrp
                     de = de * fgrp
                     dc = dc * fgrp
                  end if
c
c     form the chain rule terms for derivative expressions
c
                  de = de / r
                  dedx = de * xr
                  dedy = de * yr
                  dedz = de * zr
                  dedxc = dc * xc
                  dedyc = dc * yc
                  dedzc = dc * zc
c
c     increment the overall energy and derivative expressions
c
                  ec = ec + e
                  dec(1,i) = dec(1,i) + dedx
                  dec(2,i) = dec(2,i) + dedy
                  dec(3,i) = dec(3,i) + dedz
                  dec(1,ic) = dec(1,ic) + dedxc
                  dec(2,ic) = dec(2,ic) + dedyc
                  dec(3,ic) = dec(3,ic) + dedzc
                  dec(1,k) = dec(1,k) - dedx
                  dec(2,k) = dec(2,k) - dedy
                  dec(3,k) = dec(3,k) - dedz
                  dec(1,kc) = dec(1,kc) - dedxc
                  dec(2,kc) = dec(2,kc) - dedyc
                  dec(3,kc) = dec(3,kc) - dedzc
c
c     increment the internal virial tensor components
c
                  vxx = xr*dedx + xc*dedxc
                  vyx = yr*dedx + yc*dedxc
                  vzx = zr*dedx + zc*dedxc
                  vyy = yr*dedy + yc*dedyc
                  vzy = zr*dedy + zc*dedyc
                  vzz = zr*dedz + zc*dedzc
                  vir(1,1) = vir(1,1) + vxx
                  vir(2,1) = vir(2,1) + vyx
                  vir(3,1) = vir(3,1) + vzx
                  vir(1,2) = vir(1,2) + vyx
                  vir(2,2) = vir(2,2) + vyy
                  vir(3,2) = vir(3,2) + vzy
                  vir(1,3) = vir(1,3) + vzx
                  vir(2,3) = vir(2,3) + vzy
                  vir(3,3) = vir(3,3) + vzz
c
c     increment the total intermolecular energy
c
                  if (.not.prime .or. molcule(i).ne.molcule(k)) then
                     einter = einter + e
                  end if
               end if
            end if
   20       continue
         end do
         if (repeat) then
            repeat = .false.
            start = kbx(ii) + 1
            stop = nlight
            goto 10
         end if
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      deallocate (xsort)
      deallocate (ysort)
      deallocate (zsort)
      return
      end
c
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine echarge1c  --  charge derivs via neighbor list  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "echarge1c" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     using a pairwise neighbor list
c
c
      subroutine echarge1c
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'group.i'
      include 'inter.i'
      include 'molcul.i'
      include 'neigh.i'
      include 'shunt.i'
      include 'usage.i'
      include 'virial.i'
      integer i,j,k
      integer ii,kk,kkk
      integer in,kn,ic,kc
      real*8 e,fgrp,de,dc
      real*8 f,fi,fik
      real*8 r,r2,rb,rb2
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 xc,yc,zc
      real*8 xic,yic,zic
      real*8 dedx,dedy,dedz
      real*8 dedxc,dedyc,dedzc
      real*8 shift,taper,trans
      real*8 dtaper,dtrans
      real*8 rc,rc2,rc3,rc4
      real*8 rc5,rc6,rc7
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      real*8, allocatable :: cscale(:)
      logical proceed,usei
      character*6 mode
c
c
c     zero out the charge interaction energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     set conversion factor, cutoff and switching coefficients
c
      f = electric / dielec
      mode = 'CHARGE'
      call switch (mode)
c
c     compute the charge interaction energy and first derivatives
c
      do ii = 1, nion
         i = iion(ii)
         in = jion(ii)
         ic = kion(ii)
         xic = x(ic)
         yic = y(ic)
         zic = z(ic)
         xi = x(i) - xic
         yi = y(i) - yic
         zi = z(i) - zic
         fi = f * pchg(ii)
         usei = (use(i) .or. use(ic))
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kkk = 1, nelst(ii)
            kk = elst(kkk,ii)
            k = iion(kk)
            kn = jion(kk)
            kc = kion(kk)
            proceed = .true.
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            if (proceed)  proceed = (usei .or. use(k) .or. use(kc))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xc = xic - x(kc)
               yc = yic - y(kc)
               zc = zic - z(kc)
               if (use_bounds)  call image (xc,yc,zc)
               rc2 = xc*xc + yc*yc + zc*zc
               if (rc2 .le. off2) then
                  xr = xc + xi - x(k) + x(kc)
                  yr = yc + yi - y(k) + y(kc)
                  zr = zc + zi - z(k) + z(kc)
                  r2 = xr*xr + yr*yr + zr*zr
                  r = sqrt(r2)
                  rb = r + ebuffer
                  rb2 = rb * rb
                  fik = fi * pchg(kk) * cscale(kn)
                  e = fik / rb
                  de = -fik / rb2
                  dc = 0.0d0
c
c     use shifted energy switching if near the cutoff distance
c
                  shift = fik / (0.5d0*(off+cut))
                  e = e - shift
                  if (rc2 .gt. cut2) then
                     rc = sqrt(rc2)
                     rc3 = rc2 * rc
                     rc4 = rc2 * rc2
                     rc5 = rc2 * rc3
                     rc6 = rc3 * rc3
                     rc7 = rc3 * rc4
                     taper = c5*rc5 + c4*rc4 + c3*rc3
     &                          + c2*rc2 + c1*rc + c0
                     dtaper = 5.0d0*c5*rc4 + 4.0d0*c4*rc3
     &                           + 3.0d0*c3*rc2 + 2.0d0*c2*rc + c1
                     trans = fik * (f7*rc7 + f6*rc6 + f5*rc5 + f4*rc4
     &                               + f3*rc3 + f2*rc2 + f1*rc + f0)
                     dtrans = fik * (7.0d0*f7*rc6 + 6.0d0*f6*rc5
     &                               + 5.0d0*f5*rc4 + 4.0d0*f4*rc3
     &                             + 3.0d0*f3*rc2 + 2.0d0*f2*rc + f1)
                     dc = (e*dtaper + dtrans) / rc
                     de = de * taper
                     e = e*taper + trans
                  end if
c
c     scale the interaction based on its group membership
c
                  if (use_group) then
                     e = e * fgrp
                     de = de * fgrp
                     dc = dc * fgrp
                  end if
c
c     form the chain rule terms for derivative expressions
c
                  de = de / r
                  dedx = de * xr
                  dedy = de * yr
                  dedz = de * zr
                  dedxc = dc * xc
                  dedyc = dc * yc
                  dedzc = dc * zc
c
c     increment the overall energy and derivative expressions
c
                  ec = ec + e
                  dec(1,i) = dec(1,i) + dedx
                  dec(2,i) = dec(2,i) + dedy
                  dec(3,i) = dec(3,i) + dedz
                  dec(1,ic) = dec(1,ic) + dedxc
                  dec(2,ic) = dec(2,ic) + dedyc
                  dec(3,ic) = dec(3,ic) + dedzc
                  dec(1,k) = dec(1,k) - dedx
                  dec(2,k) = dec(2,k) - dedy
                  dec(3,k) = dec(3,k) - dedz
                  dec(1,kc) = dec(1,kc) - dedxc
                  dec(2,kc) = dec(2,kc) - dedyc
                  dec(3,kc) = dec(3,kc) - dedzc
c
c     increment the internal virial tensor components
c
                  vxx = xr*dedx + xc*dedxc
                  vyx = yr*dedx + yc*dedxc
                  vzx = zr*dedx + zc*dedxc
                  vyy = yr*dedy + yc*dedyc
                  vzy = zr*dedy + zc*dedyc
                  vzz = zr*dedz + zc*dedzc
                  vir(1,1) = vir(1,1) + vxx
                  vir(2,1) = vir(2,1) + vyx
                  vir(3,1) = vir(3,1) + vzx
                  vir(1,2) = vir(1,2) + vyx
                  vir(2,2) = vir(2,2) + vyy
                  vir(3,2) = vir(3,2) + vzy
                  vir(1,3) = vir(1,3) + vzx
                  vir(2,3) = vir(2,3) + vzy
                  vir(3,3) = vir(3,3) + vzz
c
c     increment the total intermolecular energy
c
                  if (molcule(i) .ne. molcule(k)) then
                     einter = einter + e
                  end if
               end if
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      return
      end
c
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine echarge1d  --  double loop Ewald charge derivs  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "echarge1d" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     using a particle mesh Ewald summation
c
c
      subroutine echarge1d
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'boxes.i'
      include 'cell.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'ewald.i'
      include 'group.i'
      include 'inter.i'
      include 'math.i'
      include 'molcul.i'
      include 'shunt.i'
      include 'usage.i'
      include 'virial.i'
      integer i,j,k
      integer ii,kk
      integer in,kn
      real*8 e,de,efix
      real*8 eintra
      real*8 f,fi,fik,fs
      real*8 r,r2,rew
      real*8 rb,rb2
      real*8 fgrp,term
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 xd,yd,zd
      real*8 erfc,erfterm
      real*8 scale,scaleterm
      real*8 dedx,dedy,dedz
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      real*8, allocatable :: cscale(:)
      logical proceed,usei
      character*6 mode
      external erfc
c
c
c     zero out the Ewald summation energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     zero out the intramolecular portion of the Ewald energy
c
      eintra = 0.0d0
c
c     set conversion factor, cutoff and switching coefficients
c
      f = electric / dielec
      mode = 'EWALD'
      call switch (mode)
c
c     compute the Ewald self-energy term over all the atoms
c
      fs = -f * aewald / sqrtpi
      do ii = 1, nion
         e = fs * pchg(ii)**2
         ec = ec + e
      end do
c
c     compute the cell dipole boundary correction term
c
      if (boundary .eq. 'VACUUM') then
         xd = 0.0d0
         yd = 0.0d0
         zd = 0.0d0
         do ii = 1, nion
            i = iion(ii)
            xd = xd + pchg(ii)*x(i)
            yd = yd + pchg(ii)*y(i)
            zd = zd + pchg(ii)*z(i)
         end do
         term = (2.0d0/3.0d0) * f * (pi/volbox)
         e = term * (xd*xd+yd*yd+zd*zd)
         ec = ec + e
         do ii = 1, nion
            i = iion(ii)
            de = 2.0d0 * term * pchg(ii)
            dedx = de * xd
            dedy = de * yd
            dedz = de * zd
            dec(1,i) = dec(1,i) + dedx
            dec(2,i) = dec(2,i) + dedy
            dec(3,i) = dec(3,i) + dedz
         end do
      end if
c
c     compute reciprocal space Ewald energy and first derivatives
c
      call ecrecip1
c
c     compute the real space Ewald energy and first derivatives
c
      do ii = 1, nion-1
         i = iion(ii)
         in = jion(ii)
         usei = use(i)
         xi = x(i)
         yi = y(i)
         zi = z(i)
         fi = f * pchg(ii)
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kk = ii+1, nion
            k = iion(kk)
            kn = jion(kk)
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            proceed = .true.
            if (proceed)  proceed = (usei .or. use(k))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xr = xi - x(k)
               yr = yi - y(k)
               zr = zi - z(k)
c
c     find energy for interactions within real space cutoff
c
               call image (xr,yr,zr)
               r2 = xr*xr + yr*yr + zr*zr
               if (r2 .le. off2) then
                  r = sqrt(r2)
                  rb = r + ebuffer
                  rb2 = rb * rb
                  fik = fi * pchg(kk)
                  rew = aewald * r
                  erfterm = erfc (rew)
                  scale = cscale(kn)
                  if (use_group)  scale = scale * fgrp
                  scaleterm = scale - 1.0d0
                  e = (fik/rb) * (erfterm+scaleterm)
                  de = -fik * ((erfterm+scaleterm)/rb2
     &                    + (2.0d0*aewald/sqrtpi)*exp(-rew**2)/r)
c
c     form the chain rule terms for derivative expressions
c
                  de = de / r
                  dedx = de * xr
                  dedy = de * yr
                  dedz = de * zr
c
c     increment the overall energy and derivative expressions
c
                  ec = ec + e
                  dec(1,i) = dec(1,i) + dedx
                  dec(2,i) = dec(2,i) + dedy
                  dec(3,i) = dec(3,i) + dedz
                  dec(1,k) = dec(1,k) - dedx
                  dec(2,k) = dec(2,k) - dedy
                  dec(3,k) = dec(3,k) - dedz
c
c     increment the internal virial tensor components
c
                  vxx = xr * dedx
                  vyx = yr * dedx
                  vzx = zr * dedx
                  vyy = yr * dedy
                  vzy = zr * dedy
                  vzz = zr * dedz
                  vir(1,1) = vir(1,1) + vxx
                  vir(2,1) = vir(2,1) + vyx
                  vir(3,1) = vir(3,1) + vzx
                  vir(1,2) = vir(1,2) + vyx
                  vir(2,2) = vir(2,2) + vyy
                  vir(3,2) = vir(3,2) + vzy
                  vir(1,3) = vir(1,3) + vzx
                  vir(2,3) = vir(2,3) + vzy
                  vir(3,3) = vir(3,3) + vzz
c
c     increment the total intramolecular energy
c
                  if (molcule(i) .eq. molcule(k)) then
                     efix = (fik/rb) * scale
                     eintra = eintra + efix
                  end if
               end if
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     for periodic boundary conditions with large cutoffs
c     neighbors must be found by the replicates method
c
      if (.not. use_replica)  return
c
c     calculate real space portion involving other unit cells
c
      do ii = 1, nion
         i = iion(ii)
         in = jion(ii)
         usei = use(i)
         xi = x(i)
         yi = y(i)
         zi = z(i)
         fi = f * pchg(ii)
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kk = ii, nion
            k = iion(kk)
            kn = jion(kk)
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            proceed = .true.
            if (proceed)  proceed = (usei .or. use(k))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               do j = 1, ncell
                  xr = xi - x(k)
                  yr = yi - y(k)
                  zr = zi - z(k)
c
c     find appropriate image and check the real space cutoff
c
                  call imager (xr,yr,zr,j)
                  r2 = xr*xr + yr*yr + zr*zr
                  if (r2 .le. off2) then
                     r = sqrt(r2)
                     rb = r + ebuffer
                     rb2 = rb * rb
                     fik = fi * pchg(kk)
                     rew = aewald * r
                     erfterm = erfc (rew)
                     scale = 1.0d0
                     if (use_group)  scale = scale * fgrp
                     if (use_polymer) then
                        if (r2 .le. polycut2) then
                           scale = scale * cscale(kn)
                        end if
                     end if
                     scaleterm = scale - 1.0d0
                     e = (fik/rb) * (erfterm+scaleterm)
                     de = -fik * ((erfterm+scaleterm)/rb2
     &                       + (2.0d0*aewald/sqrtpi)*exp(-rew**2)/r)
c
c     form the chain rule terms for derivative expressions
c
                     de = de / r
                     dedx = de * xr
                     dedy = de * yr
                     dedz = de * zr
c
c     increment the energy and gradient values
c
                     if (i .eq. k)  e = 0.5d0 * e
                     ec = ec + e
                     dec(1,i) = dec(1,i) + dedx
                     dec(2,i) = dec(2,i) + dedy
                     dec(3,i) = dec(3,i) + dedz
                     if (i .ne. k) then
                        dec(1,k) = dec(1,k) - dedx
                        dec(2,k) = dec(2,k) - dedy
                        dec(3,k) = dec(3,k) - dedz
                     end if
c
c     increment the internal virial tensor components
c
                     vxx = xr * dedx
                     vyx = yr * dedx
                     vzx = zr * dedx
                     vyy = yr * dedy
                     vzy = zr * dedy
                     vzz = zr * dedz
                     vir(1,1) = vir(1,1) + vxx
                     vir(2,1) = vir(2,1) + vyx
                     vir(3,1) = vir(3,1) + vzx
                     vir(1,2) = vir(1,2) + vyx
                     vir(2,2) = vir(2,2) + vyy
                     vir(3,2) = vir(3,2) + vzy
                     vir(1,3) = vir(1,3) + vzx
                     vir(2,3) = vir(2,3) + vzy
                     vir(3,3) = vir(3,3) + vzz
                  end if
               end do
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     intermolecular energy is total minus intramolecular part
c
      einter = einter + ec - eintra
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine echarge1e  --  Ewald charge derivs via lights  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "echarge1e" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     using a particle mesh Ewald summation and the method of lights
c
c
      subroutine echarge1e
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'boxes.i'
      include 'cell.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'ewald.i'
      include 'group.i'
      include 'inter.i'
      include 'light.i'
      include 'math.i'
      include 'molcul.i'
      include 'shunt.i'
      include 'usage.i'
      include 'virial.i'
      integer i,j,k
      integer ii,kk,in,kn
      integer kgy,kgz,kmap
      integer start,stop
      real*8 e,de,efix
      real*8 eintra
      real*8 f,fi,fik,fs
      real*8 r,r2,rew
      real*8 rb,rb2
      real*8 fgrp,term
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 xd,yd,zd
      real*8 erfc,erfterm
      real*8 scale,scaleterm
      real*8 dedx,dedy,dedz
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      real*8, allocatable :: cscale(:)
      real*8, allocatable :: xsort(:)
      real*8, allocatable :: ysort(:)
      real*8, allocatable :: zsort(:)
      logical proceed,usei
      logical prime,repeat
      character*6 mode
      external erfc
c
c
c     zero out the Ewald summation energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
      allocate (xsort(8*n))
      allocate (ysort(8*n))
      allocate (zsort(8*n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     zero out the intramolecular portion of the Ewald energy
c
      eintra = 0.0d0
c
c     set conversion factor, cutoff and switching coefficients
c
      f = electric / dielec
      mode = 'EWALD'
      call switch (mode)
c
c     compute the Ewald self-energy term over all the atoms
c
      fs = -f * aewald / sqrtpi
      do ii = 1, nion
         e = fs * pchg(ii)**2
         ec = ec + e
      end do
c
c     compute the cell dipole boundary correction term
c
      if (boundary .eq. 'VACUUM') then
         xd = 0.0d0
         yd = 0.0d0
         zd = 0.0d0
         do ii = 1, nion
            i = iion(ii)
            xd = xd + pchg(ii)*x(i)
            yd = yd + pchg(ii)*y(i)
            zd = zd + pchg(ii)*z(i)
         end do
         term = (2.0d0/3.0d0) * f * (pi/volbox)
         e = term * (xd*xd+yd*yd+zd*zd)
         ec = ec + e
         do ii = 1, nion
            i = iion(ii)
            de = 2.0d0 * term * pchg(ii)
            dedx = de * xd
            dedy = de * yd
            dedz = de * zd
            dec(1,i) = dec(1,i) + dedx
            dec(2,i) = dec(2,i) + dedy
            dec(3,i) = dec(3,i) + dedz
         end do
      end if
c
c     compute reciprocal space Ewald energy and first derivatives
c
      call ecrecip1
c
c     compute the real space portion of the Ewald summation;
c     transfer the interaction site coordinates to sorting arrays
c
      do i = 1, nion
         k = iion(i)
         xsort(i) = x(k)
         ysort(i) = y(k)
         zsort(i) = z(k)
      end do
c
c     use the method of lights to generate neighbors
c
      call lights (off,nion,xsort,ysort,zsort)
c
c     loop over all atoms computing the interactions
c
      do ii = 1, nion
         i = iion(ii)
         in = jion(ii)
         xi = xsort(rgx(ii))
         yi = ysort(rgy(ii))
         zi = zsort(rgz(ii))
         fi = f * pchg(ii)
         usei = (use(i))
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     loop over method of lights neighbors of current atom
c
         if (kbx(ii) .le. kex(ii)) then
            repeat = .false.
            start = kbx(ii) + 1
            stop = kex(ii)
         else
            repeat = .true.
            start = 1
            stop = kex(ii)
         end if
   10    continue
         do j = start, stop
            kk = locx(j)
            kgy = rgy(kk)
            if (kby(ii) .le. key(ii)) then
               if (kgy.lt.kby(ii) .or. kgy.gt.key(ii))  goto 20
            else
               if (kgy.lt.kby(ii) .and. kgy.gt.key(ii))  goto 20
            end if
            kgz = rgz(kk)
            if (kbz(ii) .le. kez(ii)) then
               if (kgz.lt.kbz(ii) .or. kgz.gt.kez(ii))  goto 20
            else
               if (kgz.lt.kbz(ii) .and. kgz.gt.kez(ii))  goto 20
            end if
            kmap = kk - ((kk-1)/nion)*nion
            k = iion(kmap)
            kn = jion(kmap)
            prime = (kk .le. nion)
c
c     decide whether to compute the current interaction
c
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            proceed = .true.
            if (proceed)  proceed = (usei .or. use(k))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xr = xi - xsort(j)
               yr = yi - ysort(kgy)
               zr = zi - zsort(kgz)
c
c     find energy for interactions within real space cutoff
c
               if (use_bounds) then
                  if (abs(xr) .gt. xcell2)  xr = xr - sign(xcell,xr)
                  if (abs(yr) .gt. ycell2)  yr = yr - sign(ycell,yr)
                  if (abs(zr) .gt. zcell2)  zr = zr - sign(zcell,zr)
                  if (monoclinic) then
                     xr = xr + zr*beta_cos
                     zr = zr * beta_sin
                  else if (triclinic) then
                     xr = xr + yr*gamma_cos + zr*beta_cos
                     yr = yr*gamma_sin + zr*beta_term
                     zr = zr * gamma_term
                  end if
               end if
               r2 = xr*xr + yr*yr + zr*zr
               if (r2 .le. off2) then
                  r = sqrt(r2)
                  rb = r + ebuffer
                  rb2 = rb * rb
                  rew = aewald * r
                  erfterm = erfc (rew)
                  scale = 1.0d0
                  if (prime)  scale = cscale(kn)
                  if (use_group)  scale = scale * fgrp
                  fik = fi * pchg(kmap)
                  if (use_polymer) then
                     if (r2 .gt. polycut2)  fik = fi * pchg(kmap)
                  end if
                  scaleterm = scale - 1.0d0
                  e = (fik/rb) * (erfterm+scaleterm)
                  de = -fik * ((erfterm+scaleterm)/rb2
     &                    + (2.0d0*aewald/sqrtpi)*exp(-rew**2)/r)
c
c     form the chain rule terms for derivative expressions
c
                  de = de / r
                  dedx = de * xr
                  dedy = de * yr
                  dedz = de * zr
c
c     increment the overall energy and derivative expressions
c
                  ec = ec + e
                  dec(1,i) = dec(1,i) + dedx
                  dec(2,i) = dec(2,i) + dedy
                  dec(3,i) = dec(3,i) + dedz
                  dec(1,k) = dec(1,k) - dedx
                  dec(2,k) = dec(2,k) - dedy
                  dec(3,k) = dec(3,k) - dedz
c
c     increment the internal virial tensor components
c
                  vxx = xr * dedx
                  vyx = yr * dedx
                  vzx = zr * dedx
                  vyy = yr * dedy
                  vzy = zr * dedy
                  vzz = zr * dedz
                  vir(1,1) = vir(1,1) + vxx
                  vir(2,1) = vir(2,1) + vyx
                  vir(3,1) = vir(3,1) + vzx
                  vir(1,2) = vir(1,2) + vyx
                  vir(2,2) = vir(2,2) + vyy
                  vir(3,2) = vir(3,2) + vzy
                  vir(1,3) = vir(1,3) + vzx
                  vir(2,3) = vir(2,3) + vzy
                  vir(3,3) = vir(3,3) + vzz
c
c     increment the total intramolecular energy
c
                  if (prime .and. molcule(i).eq.molcule(k)) then
                     efix = (fik/rb) * scale
                     eintra = eintra + efix
                  end if
               end if
            end if
   20       continue
         end do
         if (repeat) then
            repeat = .false.
            start = kbx(ii) + 1
            stop = nlight
            goto 10
         end if
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     intermolecular energy is total minus intramolecular part
c
      einter = einter + ec - eintra
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      deallocate (xsort)
      deallocate (ysort)
      deallocate (zsort)
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine echarge1f  --  Ewald charge derivs via list  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "echarge1f" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     using a particle mesh Ewald summation and a neighbor list
c
c
      subroutine echarge1f
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'boxes.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'ewald.i'
      include 'group.i'
      include 'inter.i'
      include 'math.i'
      include 'molcul.i'
      include 'neigh.i'
      include 'shunt.i'
      include 'usage.i'
      include 'virial.i'
      integer i,j,k
      integer ii,kk,kkk
      integer in,kn
      real*8 e,de,efix
      real*8 eintra
      real*8 f,fi,fik,fs
      real*8 r,r2,rew
      real*8 rb,rb2
      real*8 fgrp,term
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 xd,yd,zd
      real*8 erfc,erfterm
      real*8 scale,scaleterm
      real*8 dedx,dedy,dedz
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      real*8 ect,eintrat
      real*8 virt(3,3)
      real*8, allocatable :: cscale(:)
      real*8, allocatable :: dect1(:,:)
      real*8, allocatable :: dect2(:,:)
      logical proceed,usei
      character*6 mode
      external erfc
c
c
c     zero out the Ewald summation energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
      allocate (dect1(3,n))
      allocate (dect2(3,n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     zero out the intramolecular portion of the Ewald energy
c
      eintra = 0.0d0
c
c     set conversion factor, cutoff and switching coefficients
c
      f = electric / dielec
      mode = 'EWALD'
      call switch (mode)
c
c     compute the Ewald self-energy term over all the atoms
c
      fs = -f * aewald / sqrtpi
      do ii = 1, nion
         e = fs * pchg(ii)**2
         ec = ec + e
      end do
c
c     compute the cell dipole boundary correction term
c
      if (boundary .eq. 'VACUUM') then
         xd = 0.0d0
         yd = 0.0d0
         zd = 0.0d0
         do ii = 1, nion
            i = iion(ii)
            xd = xd + pchg(ii)*x(i)
            yd = yd + pchg(ii)*y(i)
            zd = zd + pchg(ii)*z(i)
         end do
         term = (2.0d0/3.0d0) * f * (pi/volbox)
         e = term * (xd*xd+yd*yd+zd*zd)
         ec = ec + e
         do ii = 1, nion
            i = iion(ii)
            de = 2.0d0 * term * pchg(ii)
            dedx = de * xd
            dedy = de * yd
            dedz = de * zd
            dec(1,i) = dec(1,i) + dedx
            dec(2,i) = dec(2,i) + dedy
            dec(3,i) = dec(3,i) + dedz
         end do
      end if
c
c     compute reciprocal space Ewald energy and first derivatives
c
      call ecrecip1
c
c     initialize local variables for OpenMP calculation
c
      ect = ec
      eintrat = eintra
      do i = 1, n
         dect1(1,i) = 0.0d0
         dect1(2,i) = 0.0d0
         dect1(3,i) = 0.0d0
         dect2(1,i) = 0.0d0
         dect2(2,i) = 0.0d0
         dect2(3,i) = 0.0d0
      end do
      do i = 1, 3
         virt(1,i) = vir(1,i)
         virt(2,i) = vir(2,i)
         virt(3,i) = vir(3,i)
      end do
c
c     set OpenMP directives for the major loop structure
c
!$OMP PARALLEL default(private) shared(nion,iion,jion,use,
!$OMP& x,y,z,f,pchg,nelst,elst,n12,n13,n14,n15,i12,i13,i14,
!$OMP& i15,c2scale,c3scale,c4scale,c5scale,use_group,off2,
!$OMP& aewald,molcule,ebuffer) firstprivate(cscale)
!$OMP& shared(ect,eintrat,dect1,dect2,virt)
!$OMP DO reduction(+:ect,eintrat,dect1,dect2,virt)
!$OMP& schedule(dynamic)
c
c     compute the real space Ewald energy and first derivatives
c
      do ii = 1, nion
         i = iion(ii)
         in = jion(ii)
         usei = use(i)
         xi = x(i)
         yi = y(i)
         zi = z(i)
         fi = f * pchg(ii)
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kkk = 1, nelst(ii)
            kk = elst(kkk,ii)
            k = iion(kk)
            kn = jion(kk)
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            proceed = .true.
            if (proceed)  proceed = (usei .or. use(k))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xr = xi - x(k)
               yr = yi - y(k)
               zr = zi - z(k)
c
c     find energy for interactions within real space cutoff
c
               call image (xr,yr,zr)
               r2 = xr*xr + yr*yr + zr*zr
               if (r2 .le. off2) then
                  r = sqrt(r2)
                  rb = r + ebuffer
                  rb2 = rb * rb
                  fik = fi * pchg(kk)
                  rew = aewald * r
                  erfterm = erfc (rew)
                  scale = cscale(kn)
                  if (use_group)  scale = scale * fgrp
                  scaleterm = scale - 1.0d0
                  e = (fik/rb) * (erfterm+scaleterm)
                  de = -fik * ((erfterm+scaleterm)/rb2
     &                    + (2.0d0*aewald/sqrtpi)*exp(-rew**2)/r)
c
c     form the chain rule terms for derivative expressions
c
                  de = de / r
                  dedx = de * xr
                  dedy = de * yr
                  dedz = de * zr
c
c     increment the overall energy and derivative expressions
c
                 ect = ect + e
                 dect1(1,i) = dect1(1,i) + dedx
                 dect1(2,i) = dect1(2,i) + dedy
                 dect1(3,i) = dect1(3,i) + dedz
                 dect2(1,k) = dect2(1,k) - dedx
                 dect2(2,k) = dect2(2,k) - dedy
                 dect2(3,k) = dect2(3,k) - dedz
c
c     increment the internal virial tensor components
c
                  vxx = xr * dedx
                  vyx = yr * dedx
                  vzx = zr * dedx
                  vyy = yr * dedy
                  vzy = zr * dedy
                  vzz = zr * dedz
                  virt(1,1) = virt(1,1) + vxx
                  virt(2,1) = virt(2,1) + vyx
                  virt(3,1) = virt(3,1) + vzx
                  virt(1,2) = virt(1,2) + vyx
                  virt(2,2) = virt(2,2) + vyy
                  virt(3,2) = virt(3,2) + vzy
                  virt(1,3) = virt(1,3) + vzx
                  virt(2,3) = virt(2,3) + vzy
                  virt(3,3) = virt(3,3) + vzz
c
c     increment the total intramolecular energy
c
                  if (molcule(i) .eq. molcule(k)) then
                     efix = (fik/rb) * scale
                     eintrat = eintrat + efix
                  end if
               end if
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     end OpenMP directives for the major loop structure
c
!$OMP END DO
!$OMP END PARALLEL
c
c     add local copies to global variables for OpenMP calculation
c
      ec = ect
      eintra = eintrat
      do i = 1, n
         dec(1,i) = dec(1,i) + dect1(1,i) + dect2(1,i)
         dec(2,i) = dec(2,i) + dect1(2,i) + dect2(2,i)
         dec(3,i) = dec(3,i) + dect1(3,i) + dect2(3,i)
      end do
      do i = 1, 3
         vir(1,i) = virt(1,i)
         vir(2,i) = virt(2,i)
         vir(3,i) = virt(3,i)
      end do
c
c     intermolecular energy is total minus intramolecular part
c
      einter = einter + ec - eintra
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      deallocate (dect1)
      deallocate (dect2)
      return
      end
c
c
c     ##################################################################
c     ##                                                              ##
c     ##  subroutine echarge1g  --  charge derivatives for smoothing  ##
c     ##                                                              ##
c     ##################################################################
c
c
c     "echarge1g" calculates the charge-charge interaction energy
c     and first derivatives with respect to Cartesian coordinates
c     for use with potential smoothing methods
c
c
      subroutine echarge1g
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'couple.i'
      include 'deriv.i'
      include 'energi.i'
      include 'group.i'
      include 'inter.i'
      include 'math.i'
      include 'molcul.i'
      include 'usage.i'
      include 'warp.i'
      integer i,j,k
      integer ii,kk
      integer in,kn
      real*8 e,de,fgrp
      real*8 r,r2,rb,rb2
      real*8 f,fi,fik
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 dedx,dedy,dedz
      real*8 erf,erfterm
      real*8 expcut,expterm
      real*8 wterm,width
      real*8 width2,width3
      real*8, allocatable :: cscale(:)
      logical proceed,usei
      external erf
c
c
c     zero out the charge interaction energy and derivatives
c
      ec = 0.0d0
      do i = 1, n
         dec(1,i) = 0.0d0
         dec(2,i) = 0.0d0
         dec(3,i) = 0.0d0
      end do
      if (nion .eq. 0)  return
c
c     perform dynamic allocation of some local arrays
c
      allocate (cscale(n))
c
c     set array needed to scale connected atom interactions
c
      do i = 1, n
         cscale(i) = 1.0d0
      end do
c
c     set the energy units conversion factor
c
      f = electric / dielec
      expcut = -50.0d0
c
c     set the extent of smoothing to be performed
c
      width = deform * diffc
      if (use_dem) then
         if (width .gt. 0.0d0)  width = 0.5d0 / sqrt(width)
      else if (use_gda) then
         wterm = sqrt(3.0d0/(2.0d0*diffc))
      end if
      width2 = width * width
      width3 = width * width2
c
c     compute the charge interaction energy and first derivatives
c
      do ii = 1, nion-1
         i = iion(ii)
         in = jion(ii)
         usei = (use(i))
         xi = x(i)
         yi = y(i)
         zi = z(i)
         fi = f * pchg(ii)
c
c     set interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = c2scale
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = c3scale
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = c4scale
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = c5scale
         end do
c
c     decide whether to compute the current interaction
c
         do kk = ii+1, nion
            k = iion(kk)
            kn = jion(kk)
            proceed = .true.
            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
            if (proceed)  proceed = (usei .or. use(k))
c
c     compute the energy contribution for this interaction
c
            if (proceed) then
               xr = xi - x(k)
               yr = yi - y(k)
               zr = zi - z(k)
               r2 = xr*xr + yr*yr + zr*zr
               r = sqrt(r2)
               rb = r + ebuffer
               rb2 = rb * rb
               fik = fi * pchg(kk) * cscale(kn)
               e = fik / rb
               de = -fik / rb2
c
c     transform the potential function via smoothing
c
               if (use_dem) then
                  if (width .gt. 0.0d0) then
                     erfterm = erf(width*rb)
                     expterm = -rb2 * width2
                     if (expterm .gt. expcut) then
                        expterm = 2.0d0*fik*width*exp(expterm)
     &                               / (sqrtpi*rb)
                     else
                        expterm = 0.0d0
                     end if
                     e = e * erfterm
                     de = de*erfterm + expterm
                  end if
               else if (use_gda) then
                  width = m2(i) + m2(k)
                  if (width .gt. 0.0d0) then
                     width = wterm / sqrt(width)
                     erfterm = erf(width*rb)
                     expterm = -rb2 * width * width
                     if (expterm .gt. expcut) then
                        expterm = 2.0d0*fik*width*exp(expterm)
     &                               / (sqrtpi*rb)
                     else
                        expterm = 0.0d0
                     end if
                     e = e * erfterm
                     de = de*erfterm + expterm
                  end if
               else if (use_tophat) then
                  if (width .gt. rb) then
                     e = fik * (3.0d0*width2-rb2) / (2.0d0*width3)
                     de = -fik * rb / width3
                  end if
               else if (use_stophat) then
                  wterm = rb + width
                  e = fik / wterm
                  de = -fik / (wterm*wterm)
               end if
c
c     scale the interaction based on its group membership
c
               if (use_group) then
                  e = e * fgrp
                  de = de * fgrp
               end if
c
c     form the chain rule terms for derivative expressions
c
               de = de / r
               dedx = de * xr
               dedy = de * yr
               dedz = de * zr
c
c     increment the overall energy and derivative expressions
c
               ec = ec + e
               dec(1,i) = dec(1,i) + dedx
               dec(2,i) = dec(2,i) + dedy
               dec(3,i) = dec(3,i) + dedz
               dec(1,k) = dec(1,k) - dedx
               dec(2,k) = dec(2,k) - dedy
               dec(3,k) = dec(3,k) - dedz
c
c     increment the total intermolecular energy
c
               if (molcule(i) .ne. molcule(k)) then
                  einter = einter + e
               end if
            end if
         end do
c
c     reset interaction scaling coefficients for connected atoms
c
         do j = 1, n12(in)
            cscale(i12(j,in)) = 1.0d0
         end do
         do j = 1, n13(in)
            cscale(i13(j,in)) = 1.0d0
         end do
         do j = 1, n14(in)
            cscale(i14(j,in)) = 1.0d0
         end do
         do j = 1, n15(in)
            cscale(i15(j,in)) = 1.0d0
         end do
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (cscale)
      return
      end
c
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine ecrecip1  --  PME recip energy and derivatives  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "ecrecip1" evaluates the reciprocal space portion of the particle
c     mesh Ewald summation energy and gradient due to partial charges
c
c     literature reference:
c
c     U. Essmann, L. Perera, M. L Berkowitz, T. Darden, H. Lee and
c     L. G. Pedersen, "A Smooth Particle Mesh Ewald Method", Journal
c     of Chemical Physics, 103, 8577-8593 (1995)
c
c     modifications for nonperiodic systems suggested by Tom Darden
c     during May 2007
c
c
      subroutine ecrecip1
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'bound.i'
      include 'boxes.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'deriv.i'
      include 'energi.i'
      include 'ewald.i'
      include 'math.i'
      include 'pme.i'
      include 'virial.i'
      integer i,j,k
      integer isite,iatm
      integer i0,j0,k0
      integer it1,it2,it3
      integer k1,k2,k3
      integer m1,m2,m3
      integer nf1,nf2,nf3
      integer nff,npoint
      integer igrd0,jgrd0,kgrd0
      real*8 e,term,expterm
      real*8 vterm,pterm
      real*8 volterm
      real*8 f,fi,denom
      real*8 hsq,struc2
      real*8 de1,de2,de3
      real*8 dn1,dn2,dn3
      real*8 t1,t2,t3
      real*8 dt1,dt2,dt3
      real*8 h1,h2,h3
      real*8 r1,r2,r3
c
c
c     zero out the particle mesh Ewald charge grid
c
      do k = 1, nfft3
         do j = 1, nfft2
            do i = 1, nfft1
               qgrid(1,i,j,k) = 0.0d0
               qgrid(2,i,j,k) = 0.0d0
            end do
         end do
      end do
c
c     get B-spline coefficients and put charges onto grid
c
      call bspline_fill
      call table_fill
      call grid_pchg
c
c     perform the 3-D FFT forward transformation
c
      call fftfront
c
c     use scalar sum to get reciprocal space energy and virial
c
      f = 0.5d0 * electric / dielec
      npoint = nfft1 * nfft2 * nfft3
      pterm = (pi/aewald)**2
      volterm = pi * volbox
      nff = nfft1 * nfft2
      nf1 = (nfft1+1) / 2
      nf2 = (nfft2+1) / 2
      nf3 = (nfft3+1) / 2
      do i = 1, npoint-1
         k3 = i/nff + 1
         j = i - (k3-1)*nff
         k2 = j/nfft1 + 1
         k1 = j - (k2-1)*nfft1 + 1
         m1 = k1 - 1
         m2 = k2 - 1
         m3 = k3 - 1
         if (k1 .gt. nf1)  m1 = m1 - nfft1
         if (k2 .gt. nf2)  m2 = m2 - nfft2
         if (k3 .gt. nf3)  m3 = m3 - nfft3
         r1 = dble(m1)
         r2 = dble(m2)
         r3 = dble(m3)
         h1 = recip(1,1)*r1 + recip(1,2)*r2 + recip(1,3)*r3
         h2 = recip(2,1)*r1 + recip(2,2)*r2 + recip(2,3)*r3
         h3 = recip(3,1)*r1 + recip(3,2)*r2 + recip(3,3)*r3
         hsq = h1*h1 + h2*h2 + h3*h3
         term = -pterm * hsq
         expterm = 0.0d0
         if (term .gt. -50.0d0) then
            denom = volterm*hsq*bsmod1(k1)*bsmod2(k2)*bsmod3(k3)
            expterm = exp(term) / denom
            if (.not. use_bounds) then
               expterm = expterm * (1.0d0-cos(pi*xbox*sqrt(hsq)))
            else if (octahedron) then
               if (mod(m1+m2+m3,2) .ne. 0)  expterm = 0.0d0
            end if
            struc2 = qgrid(1,k1,k2,k3)**2 + qgrid(2,k1,k2,k3)**2
            e = f * expterm * struc2
            ec = ec + e
            vterm = (2.0d0/hsq) * (1.0d0-term) * e
            vir(1,1) = vir(1,1) + h1*h1*vterm - e
            vir(2,1) = vir(2,1) + h1*h2*vterm
            vir(3,1) = vir(3,1) + h1*h3*vterm
            vir(1,2) = vir(1,2) + h2*h1*vterm
            vir(2,2) = vir(2,2) + h2*h2*vterm - e
            vir(3,2) = vir(3,2) + h2*h3*vterm
            vir(1,3) = vir(1,3) + h3*h1*vterm
            vir(2,3) = vir(2,3) + h3*h2*vterm
            vir(3,3) = vir(3,3) + h3*h3*vterm - e
         end if
         qgrid(1,k1,k2,k3) = expterm * qgrid(1,k1,k2,k3)
         qgrid(2,k1,k2,k3) = expterm * qgrid(2,k1,k2,k3)
      end do
c
c     account for the zeroth grid point for a finite system
c
      if (.not. use_bounds) then
         expterm = 0.5d0 * pi / xbox
         struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2
         e = f * expterm * struc2
         ec = ec + e
         qgrid(1,1,1,1) = expterm * qgrid(1,1,1,1)
         qgrid(2,1,1,1) = expterm * qgrid(2,1,1,1)
      end if
c
c     perform the 3-D FFT backward transformation
c
      call fftback
c
c     get first derivatives of the reciprocal space energy
c
      f = electric / dielec
      dn1 = dble(nfft1)
      dn2 = dble(nfft2)
      dn3 = dble(nfft3)
      do isite = 1, nion
         iatm = iion(isite)
         igrd0 = igrid(1,iatm)
         jgrd0 = igrid(2,iatm)
         kgrd0 = igrid(3,iatm)
         fi = f * pchg(isite)
         de1 = 0.0d0
         de2 = 0.0d0
         de3 = 0.0d0
         k0 = kgrd0
         do it3 = 1, bsorder
            k0 = k0 + 1
            k = k0 + 1 + (nfft3-sign(nfft3,k0))/2
            t3 = thetai3(1,it3,iatm)
            dt3 = dn3 * thetai3(2,it3,iatm)
            j0 = jgrd0
            do it2 = 1, bsorder
               j0 = j0 + 1
               j = j0 + 1 + (nfft2-sign(nfft2,j0))/2
               t2 = thetai2(1,it2,iatm)
               dt2 = dn2 * thetai2(2,it2,iatm)
               i0 = igrd0
               do it1 = 1, bsorder
                  i0 = i0 + 1
                  i = i0 + 1 + (nfft1-sign(nfft1,i0))/2
                  t1 = thetai1(1,it1,iatm)
                  dt1 = dn1 * thetai1(2,it1,iatm)
                  term = qgrid(1,i,j,k)
                  de1 = de1 + term*dt1*t2*t3
                  de2 = de2 + term*dt2*t1*t3
                  de3 = de3 + term*dt3*t1*t2
               end do
            end do
         end do
         dec(1,iatm) = dec(1,iatm) + fi*(recip(1,1)*de1+recip(1,2)*de2
     &                                      +recip(1,3)*de3)
         dec(2,iatm) = dec(2,iatm) + fi*(recip(2,1)*de1+recip(2,2)*de2
     &                                      +recip(2,3)*de3)
         dec(3,iatm) = dec(3,iatm) + fi*(recip(3,1)*de1+recip(3,2)*de2
     &                                      +recip(3,3)*de3)
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


module l2m_test
    use fruit
    use global_defs
    implicit none

contains

    subroutine test_l2m_getlvl2minors
        integer, parameter :: n = 4
        integer, parameter :: nv1 = 0
        integer, parameter :: nv2 = 0
        real(dp) :: csc(n+nv1, n+nv2)
        real(dp) :: l2minor(n**4)
        real(dp), parameter :: tol = 1.E-12_dp
        integer :: i
        external :: getlvl2minors_lu

        csc = reshape([ (i, i=1, (n+nv1)*(n+nv2)) ], shape(csc), order=[2, 1])

        call getlvl2minors_lu(n, nv2, csc, n+nv1, l2minor, n**3)
        call assert_equals(-4.0_dp, l2minor(18), tol, 'Level 2 minors')
        call assert_equals(-4.0_dp, l2minor(79), tol, 'Level 2 minors')
        call assert_equals(-4.0_dp, l2minor(226), tol, 'Level 2 minors')
        call assert_equals(-4.0_dp, l2minor(191), tol, 'Level 2 minors')
    end subroutine test_l2m_getlvl2minors

end module l2m_test


text:
Debras Electrolysis offers body and nostril waxing, facials and Eyelash extensions by Novalash as well as permanent hair removal. Also a certified trainer for Novalash training and certifying beauty and medical professionals in the art of applying eyelash extensions.

Offering:
  Eyeliner
  Lash Enhancement
  Eye Brows
  Lip Liner
  Permanent Hair Removal

  To learn about other businesses in Davis that provide professional aesthetic services, visit our Beauty Salons and Barber Shops page.

20081012 12:55:59 nbsp Debras Electrolysis also offers body and nostril waxing, facials and Eyelash extensions by Novalash as well as permanent hair removal. Also a certified trainer for Novalash training and certifying beauty and medical professionals in the art of applying eyelash extensions.    Users/Debrakontos

20090413 17:51:31 nbsp Does anybody know if she accepts debit/credit card? Users/JAILEDAN
  Yes. My wife goes to her every week and says she accepts both debit and credit. Users/condemned2bfree

20090522 14:47:04 nbsp I go here regularly for electrolysis, so far the results have been great. It helps that Deborah is such a sweet person and easy to talk to, this makes the electrolysis process seem much less painful. Users/Zeeba

20090714 21:43:41 nbsp Thanks for the reply, condemned2bfree.. sorry took so long to comment back!

I definitely recommend Deborah for a Brazilian wax!  It was my first time when I had my wax done by her, and she just made it soooooo comfortable and easy.  She used hard wax that did not need waxing strips.  Once it hardened, she peeled it off without pulling the skin.  It only stuck to the hair!!  Deborah also explained the process really well, and she was a very interesting person to talk to.  I feel bad that I havent returned since my first and last visit, because... well.. its kinda expensive for me about $70 for a full monty (WELL WORTH IT, I GUARANTEE!!!).  It depends on the volume of hair you want to get rid of.  :)  I shall see her again soon!   Users/JAILEDAN


text:
	PROGRAM MAKEHOURLY
C
C	TO COMBINE TWO FILES WITH DIFFERENT TIME SCALES
C
	REAL*8 SCETT
	character*79 junk
	INTEGER*4 YYYY,HR
C
	OPEN(UNIT=44,FILE='DATA_A:[KELLOGG]WI_MFI_199906.1HR',STATUS='OLD',
     1		READONLY)
 1001	FORMAT(A)
	DO N = 1,97
C	DO N = 1,98
	  READ(44,1001) JUNK
	ENDDO
C	PRINT 1001, JUNK
C
	OPEN(UNIT=46,FILE='DATA_A:[KELLOGG]WI_3DP_199906.DAT',STATUS='OLD',
     1		READONLY)
	DO N = 1,74
	  READ(46,1001) JUNK
	ENDDO
	PRINT 1001, JUNK
	if(1) stop
C
 100	CONTINUE
	READ(44,1077,ERR=100,END=300) MDAY,MO,YYYY,HR,MIN,SEC,BMAG,RRE
 1077	  FORMAT(1X,I2,1X,I2,1X,I4,I3,1X,I2,1X,F6.3,3F17.5)
	  SCETT = 3600.D00*HR + 60.D00*MIN + SEC
	PRINT 1077,MDAY,MO,YYYY,HR,MIN,SEC,BMAG,RRE
	PRINT*,'44 SEC OF DAY',SCETT
	READ(46,1077,ERR=100,END=300) MDAY,MO,YYYY,HR,MIN,SEC,TE,DENS
     1		,TI,VX,VY,VZ
	  SCETT = 3600.D00*HR + 60.D00*MIN + SEC
	PRINT 1077,MDAY,MO,YYYY,HR,MIN,SEC,TE,DENS,TI,VX,VY,VZ
	PRINT*,'46 SEC OF DAY',SCETT
 300	CONTINUE
	STOP
	END

ext: for


text:
! Example of an approximate stencil specification
!
! In situations where the specification language cannot exactly
! describe the pattern an approximate specification can be given as an
! upper and lower bound using the atMost and atLeast modifiers
program approx
implicit none

integer i, n
parameter (n = 3)
real a(0:n)

do i = 0, n
 != stencil readOnce, atLeast, pointed(dim=1) :: a
 != stencil readOnce, atMost, forward(depth=4, dim=1) :: a
 a(i) = a(i) + a(i+4)
end do

end program

max_stars_repo_name: raehik/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: raehik/camfort


text:
  subroutine constructGeneralData( gsmObj )

! ==============================================================================
!
! This function constructs the general physics data used by GSM in the generalPhysicsData
! data type
!
!
! Written by CMJ, XCP-3, 01/2019
!
! ==============================================================================

    use molnixClass, only: newMolnix
    use fissionBarrierClass, only: newFissionBarrier
    use evaporationDataClass, only: newEvaporationData

    implicit none
    class(GSM), intent(inout) :: gsmObj

! ==============================================================================

    gsmObj%genData%numConstructionErrors = 0

    ! Molnix object:
    gsmObj%genData%molnixEnergies = newMolnix( &
         & gsmObj%options%molnixOpts, &
         & gsmObj%io%print )
    if ( .not.gsmObj%genData%molnixEnergies%properlyConstructed() ) then
       write(gsmObj%io%message, 1000) "The Molnix object failed to properly construct."
       call gsmObj%io%print(1, 1, gsmObj%io%message)
       gsmObj%genData%numConstructionErrors = gsmObj%genData%numConstructionErrors + 1
    end if
    gsmObj%options%molnixOpts = gsmObj%genData%molnixEnergies%queryOptions()

    ! Fission Barrier object:
    gsmObj%genData%fissBarrier = newFissionBarrier( &
         & gsmObj%genData%molnixEnergies, &
         & gsmObj%options%fissBarOpts, gsmObj%io%print )
    if ( .not.gsmObj%genData%fissBarrier%properlyConstructed() ) then
       write(gsmObj%io%message, 1000) "The Fission Barrier object failed to properly construct."
       call gsmObj%io%print(1, 1, gsmObj%io%message)
       gsmObj%genData%numConstructionErrors = gsmObj%genData%numConstructionErrors + 1
    end if
    gsmObj%options%fissBarOpts = gsmObj%genData%fissBarrier%queryOptions()

    ! Evaporation Data object:
    gsmObj%genData%evap = newEvaporationData( &
         & gsmObj%options%evapDataOpts, &
         & gsmObj%io%print )
    if ( .not.gsmObj%genData%evap%properlyConstructed() ) then
       write(gsmObj%io%message, 1000) "The EvaporationData object failed to properly construct."
       call gsmObj%io%print(1, 1, gsmObj%io%message)
       gsmObj%genData%numConstructionErrors = gsmObj%genData%numConstructionErrors + 1
    end if
    gsmObj%options%evapDataOpts = gsmObj%genData%evap%queryOptions()

    ! Check how many errors occurred:
    if ( gsmObj%genData%numConstructionErrors /= 0 ) then
       write(gsmObj%io%message, 1100) gsmObj%genData%numConstructionErrors
       call gsmObj%io%print(1, 1, gsmObj%io%message)
    end if

    return
! ==============================================================================
1000 format(A)
1100 format(i2, "  data objects failed to construct within GSM.")
! ==============================================================================
  end subroutine constructGeneralData

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense