ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1707 by gezelter, Thu Nov 4 16:20:37 2004 UTC vs.
Revision 2357 by chuckv, Wed Oct 12 20:18:17 2005 UTC

# Line 1 | Line 1
1 + !!
2 + !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + !!
4 + !! The University of Notre Dame grants you ("Licensee") a
5 + !! non-exclusive, royalty free, license to use, modify and
6 + !! redistribute this software in source and binary code form, provided
7 + !! that the following conditions are met:
8 + !!
9 + !! 1. Acknowledgement of the program authors must be made in any
10 + !!    publication of scientific results based in part on use of the
11 + !!    program.  An acceptable form of acknowledgement is citation of
12 + !!    the article in which the program was described (Matthew
13 + !!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + !!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + !!    Parallel Simulation Engine for Molecular Dynamics,"
16 + !!    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + !!
18 + !! 2. Redistributions of source code must retain the above copyright
19 + !!    notice, this list of conditions and the following disclaimer.
20 + !!
21 + !! 3. Redistributions in binary form must reproduce the above copyright
22 + !!    notice, this list of conditions and the following disclaimer in the
23 + !!    documentation and/or other materials provided with the
24 + !!    distribution.
25 + !!
26 + !! This software is provided "AS IS," without a warranty of any
27 + !! kind. All express or implied conditions, representations and
28 + !! warranties, including any implied warranty of merchantability,
29 + !! fitness for a particular purpose or non-infringement, are hereby
30 + !! excluded.  The University of Notre Dame and its licensors shall not
31 + !! be liable for any damages suffered by licensee as a result of
32 + !! using, modifying or distributing the software or its
33 + !! derivatives. In no event will the University of Notre Dame or its
34 + !! licensors be liable for any lost revenue, profit or data, or for
35 + !! direct, indirect, special, consequential, incidental or punitive
36 + !! damages, however caused and regardless of the theory of liability,
37 + !! arising out of the use of or inability to use software, even if the
38 + !! University of Notre Dame has been advised of the possibility of
39 + !! such damages.
40 + !!
41 +
42 +
43   module shapes
44  
45    use force_globals
# Line 13 | Line 55 | module shapes
55    implicit none
56  
57    PRIVATE
58 <  
58 > #define __FORTRAN90
59 > #include "UseTheForce/DarkSide/fInteractionMap.h"
60    INTEGER, PARAMETER:: CHEBYSHEV_TN = 1
61    INTEGER, PARAMETER:: CHEBYSHEV_UN = 2
62    INTEGER, PARAMETER:: LAGUERRE     = 3
# Line 26 | Line 69 | module shapes
69    public :: do_shape_pair
70    public :: newShapeType
71    public :: complete_Shape_FF
72 +  public :: destroyShapeTypes
73 +  public :: getShapeCut
74  
30 –
75    type, private :: Shape
76       integer :: atid
77       integer :: nContactFuncs
# Line 51 | Line 95 | module shapes
95       real ( kind = dp )  :: epsilon
96       real ( kind = dp )  :: sigma
97    end type Shape
98 <  
98 >
99    type, private :: ShapeList
100       integer :: n_shapes = 0
101       integer :: currentShape = 0
102 <     type (Shape), pointer :: Shapes(:)      => null()
103 <     integer, pointer      :: atidToShape(:) => null()
102 >     type(Shape), pointer :: Shapes(:)      => null()
103 >     integer, pointer     :: atidToShape(:) => null()
104    end type ShapeList
105    
106    type(ShapeList), save :: ShapeMap
107 <
107 >  
108    integer :: lmax
109 <
109 >  
110   contains  
111 <
111 >  
112    subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
113         ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
114         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
115         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
116         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
117 <       myATID, status)
118 <
117 >       c_ident, status)
118 >    
119      integer :: nContactFuncs
120      integer :: nRangeFuncs
121      integer :: nStrengthFuncs
122      integer :: shape_ident
123      integer :: status
124 +    integer :: c_ident
125      integer :: myATID
126      integer :: bigL
127      integer :: bigM
# Line 102 | Line 147 | contains  
147  
148         call getMatchingElementList(atypes, "is_Shape", .true., &
149              nShapeTypes, MatchList)
150 <      
150 >
151         call getMatchingElementList(atypes, "is_LennardJones", .true., &
152              nLJTypes, MatchList)
153 <      
153 >
154         ShapeMap%n_shapes = nShapeTypes + nLJTypes
155 <      
155 >
156         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
157 <      
157 >
158         ntypes = getSize(atypes)
159 <      
160 <       allocate(ShapeMap%atidToShape(0:ntypes))
159 >
160 >       allocate(ShapeMap%atidToShape(ntypes))
161      end if
162 <    
162 >
163      ShapeMap%currentShape = ShapeMap%currentShape + 1
164      current = ShapeMap%currentShape
165  
# Line 125 | Line 170 | contains  
170         return
171      endif
172  
173 <    call getElementProperty(atypes, myATID, 'c_ident', me)
173 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
174  
175 <    ShapeMap%atidToShape(me)                         = current
176 <    ShapeMap%Shapes(current)%atid                    = me
175 >    ShapeMap%atidToShape(myATID)                     = current
176 >    ShapeMap%Shapes(current)%atid                    = myATID
177      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
178      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
179      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 147 | Line 192 | contains  
192  
193      bigL = -1
194      bigM = -1
195 <    
195 >
196      do j = 1, ShapeMap%Shapes(current)%nContactFuncs
197         if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
198            bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
# Line 185 | Line 230 | contains  
230      type(Shape), intent(inout) :: myShape
231      integer, intent(out) :: stat
232      integer :: alloc_stat
233 <
233 >
234      stat = 0
235      if (associated(myShape%contactFuncLValue)) then
236         deallocate(myShape%contactFuncLValue)
# Line 289 | Line 334 | contains  
334      return
335  
336    end subroutine allocateShape
337 <    
337 >
338    subroutine complete_Shape_FF(status)
339      integer :: status
340      integer :: i, j, l, m, lm, function_type
341      real(kind=dp) :: thisDP, sigma
342 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
342 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
343      logical :: thisProperty
344  
345      status = 0
# Line 303 | Line 348 | contains  
348         status = -1
349         return
350      end if
351 <    
351 >
352      nAtypes = getSize(atypes)
353  
354      if (nAtypes == 0) then
355         status = -1
356         return
357 <    end if
357 >    end if      
358  
359      ! atypes comes from c side
360 <    do i = 0, nAtypes
361 <      
362 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
363 <      
360 >    do i = 1, nAtypes
361 >    
362 >       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
363 >       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
364 >        
365         if (thisProperty) then
320 –          
366            ShapeMap%currentShape = ShapeMap%currentShape + 1
367            current = ShapeMap%currentShape
368  
369 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
370 <          ShapeMap%atidToShape(thisIP) = current
326 <          ShapeMap%Shapes(current)%atid = thisIP
369 >          ShapeMap%atidToShape(myATID) = current
370 >          ShapeMap%Shapes(current)%atid = myATID
371  
372            ShapeMap%Shapes(current)%isLJ = .true.
373  
374 <          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
375 <          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
376 <          
374 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
375 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
376 >
377         endif
378 <      
378 >
379      end do
380  
381      haveShapeMap = .true.
382 <    
382 >
383 > !    do i = 1, ShapeMap%n_shapes
384 > !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
385 > !    end do
386 >
387    end subroutine complete_Shape_FF
388 <    
388 >
389 >  function getShapeCut(atomID) result(cutValue)
390 >    integer, intent(in) :: atomID
391 >    real(kind=dp) :: cutValue, whoopdedoo
392 >
393 >    !! this is just a placeholder for a cutoff value, hopefully we'll
394 >    !! develop a method to calculate a sensible value
395 >    whoopdedoo = 9.0_dp
396 >
397 >    cutValue = whoopdedoo
398 >
399 >  end function getShapeCut
400 >
401    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
402         pot, A, f, t, do_pot)
403 <    
403 >
404      INTEGER, PARAMETER:: LMAX         = 64
405      INTEGER, PARAMETER:: MMAX         = 64
406  
# Line 348 | Line 408 | contains  
408      real (kind=dp), intent(inout) :: rij, r2
409      real (kind=dp), dimension(3), intent(in) :: d
410      real (kind=dp), dimension(3), intent(inout) :: fpair
411 <    real (kind=dp) :: pot, vpair, sw
411 >    real (kind=dp) :: pot, vpair, sw, dswdr
412      real (kind=dp), dimension(9,nLocal) :: A
413      real (kind=dp), dimension(3,nLocal) :: f
414      real (kind=dp), dimension(3,nLocal) :: t
# Line 359 | Line 419 | contains  
419      integer :: l, m, lm, id1, id2, localError, function_type
420      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
421      real (kind=dp) :: coeff
422 +    real (kind=dp) :: pot_temp
423  
424      real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
425      real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
# Line 410 | Line 471 | contains  
471      real (kind=dp) :: dsduxi, dsduyi, dsduzi
472      real (kind=dp) :: dsdxj, dsdyj, dsdzj
473      real (kind=dp) :: dsduxj, dsduyj, dsduzj
474 <    
474 >
475      real (kind=dp) :: depsdxi, depsdyi, depsdzi
476      real (kind=dp) :: depsduxi, depsduyi, depsduzi
477      real (kind=dp) :: depsdxj, depsdyj, depsdzj
# Line 437 | Line 498 | contains  
498      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
499      real (kind=dp) :: fxradial, fyradial, fzradial
500  
501 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
502 +
503      real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
504      real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
505      real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
# Line 447 | Line 510 | contains  
510         return      
511      endif
512  
450 –    write(*,*) rij, r2, d(1), d(2), d(3)
451 –    write(*,*) 'before, atom1, 2 = ', atom1, atom2
452 –    write(*,*) 'f1 = ', f(1,atom1), f(2,atom1), f(3,atom1)
453 –    write(*,*) 'f2 = ', f(1,atom2), f(2,atom2), f(3,atom2)
454 –    write(*,*) 't1 = ', t(1,atom1), t(2,atom1), t(3,atom1)
455 –    write(*,*) 't2 = ', t(1,atom2), t(2,atom2), t(3,atom2)
456 –    
513      !! We assume that the rotation matrices have already been calculated
514      !! and placed in the A array.
459 –
515      r3 = r2*rij
516      r5 = r3*r2
517 <    
517 >
518      drdxi = -d(1) / rij
519      drdyi = -d(2) / rij
520      drdzi = -d(3) / rij
521 +    drduxi = 0.0d0
522 +    drduyi = 0.0d0
523 +    drduzi = 0.0d0
524  
525      drdxj = d(1) / rij
526      drdyj = d(2) / rij
527      drdzj = d(3) / rij
528 <    
528 >    drduxj = 0.0d0
529 >    drduyj = 0.0d0
530 >    drduzj = 0.0d0
531 >
532      ! find the atom type id (atid) for each atom:
533   #ifdef IS_MPI
534      atid1 = atid_Row(atom1)
# Line 480 | Line 541 | contains  
541      ! use the atid to find the shape type (st) for each atom:
542      st1 = ShapeMap%atidToShape(atid1)
543      st2 = ShapeMap%atidToShape(atid2)
544 +    
545 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
546  
547      if (ShapeMap%Shapes(st1)%isLJ) then
548  
# Line 509 | Line 572 | contains  
572   #ifdef IS_MPI
573         ! rotate the inter-particle separation into the two different
574         ! body-fixed coordinate systems:
575 <      
575 >
576         xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
577         yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
578         zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
579 <      
579 >
580   #else
581         ! rotate the inter-particle separation into the two different
582         ! body-fixed coordinate systems:
583 <      
583 >
584         xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
585         yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
586         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
524 –      
525 – #endif
587  
588 + #endif
589 +       xihat = xi / rij
590 +       yihat = yi / rij
591 +       zihat = zi / rij
592         xi2 = xi*xi
593         yi2 = yi*yi
594         zi2 = zi*zi            
# Line 535 | Line 600 | contains  
600         dctidx = - zi * xi / r3
601         dctidy = - zi * yi / r3
602         dctidz = 1.0d0 / rij - zi2 / r3
603 <       dctidux =  yi / rij + (zi * yi) / r3
604 <       dctiduy = -xi / rij - (zi * xi) / r3
605 <       dctiduz = 0.0d0
603 >       dctidux = yi / rij ! - (zi * xi2) / r3
604 >       dctiduy = -xi / rij !- (zi * yi2) / r3
605 >       dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3
606  
607         ! this is an attempt to try to truncate the singularity when
608         ! sin(theta) is near 0.0:
# Line 547 | Line 612 | contains  
612            proji = sqrt(rij * 1.0d-12)
613            dcpidx = 1.0d0 / proji
614            dcpidy = 0.0d0
615 <          dcpidux = 0.0d0
616 <          dcpiduy = zi / proji
615 >          dcpidux = xi / proji
616 >          dcpiduy = 0.0d0
617            dspidx = 0.0d0
618            dspidy = 1.0d0 / proji
619 <          dspidux = -zi / proji
620 <          dspiduy = 0.0d0
619 >          dspidux = 0.0d0
620 >          dspiduy = yi / proji
621         else
622            proji = sqrt(xi2 + yi2)
623            proji3 = proji*proji*proji
624            dcpidx = 1.0d0 / proji - xi2 / proji3
625            dcpidy = - xi * yi / proji3
626 <          dcpidux = xi * yi * zi / proji3
627 <          dcpiduy = zi / proji - xi2 * zi / proji3
626 >          dcpidux = xi / proji - (xi2 * xi) / proji3
627 >          dcpiduy = - (xi * yi2) / proji3
628            dspidx = - xi * yi / proji3
629            dspidy = 1.0d0 / proji - yi2 / proji3
630 <          dspidux = -zi / proji + yi2 * zi / proji3
631 <          dspiduy = - xi * yi * zi / proji3
630 >          dspidux = - (yi * xi2) / proji3
631 >          dspiduy = yi / proji - (yi2 * yi) / proji3
632         endif
633 <      
633 >
634         cpi = xi / proji
635         dcpidz = 0.0d0
636 <       dcpiduz = -yi / proji
637 <      
636 >       dcpiduz = 0.0d0
637 >
638         spi = yi / proji
639         dspidz = 0.0d0
640 <       dspiduz = xi / proji
576 <       write(*,*) 'before lmloop', cpi, dcpidx, dcpidux
640 >       dspiduz = 0.0d0
641  
642         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
643              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 583 | Line 647 | contains  
647              CHEBYSHEV_TN, tm_i, dtm_i)
648         call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
649              CHEBYSHEV_UN, um_i, dum_i)
650 <      
650 >
651         sigma_i = 0.0d0
652         s_i = 0.0d0
653         eps_i = 0.0d0
# Line 617 | Line 681 | contains  
681               dPhuncdX = coeff * dtm_i(m) * dcpidx
682               dPhuncdY = coeff * dtm_i(m) * dcpidy
683               dPhuncdZ = coeff * dtm_i(m) * dcpidz
684 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
684 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
685               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
686               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
687            else
# Line 631 | Line 695 | contains  
695            endif
696  
697            sigma_i = sigma_i + plm_i(m,l)*Phunc
698 <
698 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux
699 > !!$          write(*,*) 'Phunc = ', Phunc
700            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
701                 Phunc * dlm_i(m,l) * dctidx
702            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
703                 Phunc * dlm_i(m,l) * dctidy
704            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
705                 Phunc * dlm_i(m,l) * dctidz
641 –          
706            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
707                 Phunc * dlm_i(m,l) * dctidux
708            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
709                 Phunc * dlm_i(m,l) * dctiduy
710            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
711                 Phunc * dlm_i(m,l) * dctiduz
712 <
712 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
713 > !!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
714 > !!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
715         end do
716  
717         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 653 | Line 719 | contains  
719            m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
720            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
721            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
656 –          
657 –         write(*,*) 'in lm loop a', coeff, dtm_i(m), dcpidx
722  
659 –
723            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
724               Phunc = coeff * tm_i(m)
725               dPhuncdX = coeff * dtm_i(m) * dcpidx
726               dPhuncdY = coeff * dtm_i(m) * dcpidy
727               dPhuncdZ = coeff * dtm_i(m) * dcpidz
728 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
728 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
729               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
730               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
731            else
# Line 676 | Line 739 | contains  
739            endif
740  
741            s_i = s_i + plm_i(m,l)*Phunc
679 –          
742  
681 –          write(*,*) 'in lm loop ', dsidx, plm_i(m,l), dPhuncdX, Phunc, dlm_i(m,l), dctidx
743            dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
744                 Phunc * dlm_i(m,l) * dctidx
745            dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
746                 Phunc * dlm_i(m,l) * dctidy
747            dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
748                 Phunc * dlm_i(m,l) * dctidz
749 <          
749 >
750            dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
751                 Phunc * dlm_i(m,l) * dctidux
752            dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 694 | Line 755 | contains  
755                 Phunc * dlm_i(m,l) * dctiduz      
756  
757         end do
758 <              
758 >
759         do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
760            l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
761            m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
762            coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
763            function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
764 <          
764 >
765            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
766               Phunc = coeff * tm_i(m)
767               dPhuncdX = coeff * dtm_i(m) * dcpidx
768               dPhuncdY = coeff * dtm_i(m) * dcpidy
769               dPhuncdZ = coeff * dtm_i(m) * dcpidz
770 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
770 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
771               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
772               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
773            else
# Line 720 | Line 781 | contains  
781            endif
782  
783            eps_i = eps_i + plm_i(m,l)*Phunc
784 <          
784 >
785            depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
786                 Phunc * dlm_i(m,l) * dctidx
787            depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
788                 Phunc * dlm_i(m,l) * dctidy
789            depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
790                 Phunc * dlm_i(m,l) * dctidz
791 <          
791 >
792            depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
793                 Phunc * dlm_i(m,l) * dctidux
794            depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 738 | Line 799 | contains  
799         end do
800  
801      endif
741 –      
742 –       ! now do j:
802  
803 +    ! now do j:
804 +
805      if (ShapeMap%Shapes(st2)%isLJ) then
806         sigma_j = ShapeMap%Shapes(st2)%sigma
807         s_j = ShapeMap%Shapes(st2)%sigma
# Line 764 | Line 825 | contains  
825         depsjduy = 0.0d0
826         depsjduz = 0.0d0
827      else
828 <      
828 >
829   #ifdef IS_MPI
830         ! rotate the inter-particle separation into the two different
831         ! body-fixed coordinate systems:
832         ! negative sign because this is the vector from j to i:
833 <      
833 >
834         xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
835         yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
836         zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
# Line 777 | Line 838 | contains  
838         ! rotate the inter-particle separation into the two different
839         ! body-fixed coordinate systems:
840         ! negative sign because this is the vector from j to i:
841 <      
841 >
842         xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
843         yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
844         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
845   #endif
846 <      
846 >
847 >       xjhat = xj / rij
848 >       yjhat = yj / rij
849 >       zjhat = zj / rij
850         xj2 = xj*xj
851         yj2 = yj*yj
852         zj2 = zj*zj
853         ctj = zj / rij
854 <      
854 >
855         if (ctj .gt. 1.0_dp) ctj = 1.0_dp
856         if (ctj .lt. -1.0_dp) ctj = -1.0_dp
857  
858         dctjdx = - zj * xj / r3
859         dctjdy = - zj * yj / r3
860         dctjdz = 1.0d0 / rij - zj2 / r3
861 <       dctjdux =  yj / rij + (zj * yj) / r3
862 <       dctjduy = -xj / rij - (zj * xj) / r3
863 <       dctjduz = 0.0d0
864 <      
861 >       dctjdux = yj / rij !- (zi * xj2) / r3
862 >       dctjduy = -xj / rij !- (zj * yj2) / r3
863 >       dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3
864 >
865         ! this is an attempt to try to truncate the singularity when
866         ! sin(theta) is near 0.0:
867  
# Line 806 | Line 870 | contains  
870            projj = sqrt(rij * 1.0d-12)
871            dcpjdx = 1.0d0 / projj
872            dcpjdy = 0.0d0
873 <          dcpjdux = 0.0d0
874 <          dcpjduy = zj / projj
873 >          dcpjdux = xj / projj
874 >          dcpjduy = 0.0d0
875            dspjdx = 0.0d0
876            dspjdy = 1.0d0 / projj
877 <          dspjdux = -zj / projj
878 <          dspjduy = 0.0d0
877 >          dspjdux = 0.0d0
878 >          dspjduy = yj / projj
879         else
880            projj = sqrt(xj2 + yj2)
881            projj3 = projj*projj*projj
882            dcpjdx = 1.0d0 / projj - xj2 / projj3
883            dcpjdy = - xj * yj / projj3
884 <          dcpjdux = xj * yj * zj / projj3
885 <          dcpjduy = zj / projj - xj2 * zj / projj3
884 >          dcpjdux = xj / projj - (xj2 * xj) / projj3
885 >          dcpjduy = - (xj * yj2) / projj3
886            dspjdx = - xj * yj / projj3
887            dspjdy = 1.0d0 / projj - yj2 / projj3
888 <          dspjdux = -zj / projj + yj2 * zj / projj3
889 <          dspjduy = - xj * yj * zj / projj3
888 >          dspjdux = - (yj * xj2) / projj3
889 >          dspjduy = yj / projj - (yj2 * yj) / projj3
890         endif
891  
892         cpj = xj / projj
893         dcpjdz = 0.0d0
894 <       dcpjduz = -yj / projj
895 <      
894 >       dcpjduz = 0.0d0
895 >
896         spj = yj / projj
897         dspjdz = 0.0d0
898 <       dspjduz = xj / projj
898 >       dspjduz = 0.0d0
899  
900 +
901 + !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
902 + !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
903         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
904              ShapeMap%Shapes(st2)%bigL, LMAX, &
905              plm_j, dlm_j)
906 <      
906 >
907         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
908              CHEBYSHEV_TN, tm_j, dtm_j)
909         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
910              CHEBYSHEV_UN, um_j, dum_j)
911 <      
911 >
912         sigma_j = 0.0d0
913         s_j = 0.0d0
914         eps_j = 0.0d0
# Line 875 | Line 942 | contains  
942               dPhuncdX = coeff * dtm_j(m) * dcpjdx
943               dPhuncdY = coeff * dtm_j(m) * dcpjdy
944               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
945 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
945 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
946               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
947               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
948            else
# Line 887 | Line 954 | contains  
954               dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
955               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
956            endif
957 <
957 >
958            sigma_j = sigma_j + plm_j(m,l)*Phunc
959 <          
959 >
960            dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
961                 Phunc * dlm_j(m,l) * dctjdx
962            dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
963                 Phunc * dlm_j(m,l) * dctjdy
964            dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
965                 Phunc * dlm_j(m,l) * dctjdz
966 <          
966 >
967            dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
968                 Phunc * dlm_j(m,l) * dctjdux
969            dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
# Line 917 | Line 984 | contains  
984               dPhuncdX = coeff * dtm_j(m) * dcpjdx
985               dPhuncdY = coeff * dtm_j(m) * dcpjdy
986               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
987 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
987 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
988               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
989               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
990            else
# Line 931 | Line 998 | contains  
998            endif
999  
1000            s_j = s_j + plm_j(m,l)*Phunc
1001 <          
1001 >
1002            dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
1003                 Phunc * dlm_j(m,l) * dctjdx
1004            dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
1005                 Phunc * dlm_j(m,l) * dctjdy
1006            dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
1007                 Phunc * dlm_j(m,l) * dctjdz
1008 <          
1008 >
1009            dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
1010                 Phunc * dlm_j(m,l) * dctjdux
1011            dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 972 | Line 1039 | contains  
1039               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1040            endif
1041  
1042 + !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1043 +
1044            eps_j = eps_j + plm_j(m,l)*Phunc
1045 <          
1045 >
1046            depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1047                 Phunc * dlm_j(m,l) * dctjdx
1048            depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1049                 Phunc * dlm_j(m,l) * dctjdy
1050            depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1051                 Phunc * dlm_j(m,l) * dctjdz
1052 <          
1052 >
1053            depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1054                 Phunc * dlm_j(m,l) * dctjdux
1055            depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 995 | Line 1064 | contains  
1064      ! phew, now let's assemble the potential energy:
1065  
1066      sigma = 0.5*(sigma_i + sigma_j)
1067 <
1067 > !    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1068      dsigmadxi = 0.5*dsigmaidx
1069      dsigmadyi = 0.5*dsigmaidy
1070      dsigmadzi = 0.5*dsigmaidz
# Line 1027 | Line 1096 | contains  
1096      dsduzj = 0.5*dsjduz
1097  
1098      eps = sqrt(eps_i * eps_j)
1099 <
1100 <    write(*,*) 'sigma, s, eps = ', sigma, s, eps
1101 <
1099 > !!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1100 > !!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1101 > !!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1102      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1103      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1104      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1043 | Line 1112 | contains  
1112      depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1113      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1114      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1046 –    
1047 –    rtdenom = rij-sigma+s
1115  
1116 <    write(*,*) 'rtdenom = ', rtdenom, ' sw = ', sw
1116 > !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1117 >
1118 > !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1119 > !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1120 > !!$
1121 > !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1122 >
1123 >    rtdenom = rij-sigma+s
1124      rt = s / rtdenom
1125  
1126 <    write(*,*) 'john' , dsdxi, rt, drdxi, dsigmadxi, rtdenom
1127 <    write(*,*) 'bigboot', dsduzj, rt, drduzj, dsigmaduzj, rtdenom
1126 >    drtdxi = (dsdxi - rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1127 >    drtdyi = (dsdyi - rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1128 >    drtdzi = (dsdzi - rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1129 >    drtduxi = (dsduxi - rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1130 >    drtduyi = (dsduyi - rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1131 >    drtduzi = (dsduzi - rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1132 >    drtdxj = (dsdxj - rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1133 >    drtdyj = (dsdyj - rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1134 >    drtdzj = (dsdzj - rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1135 >    drtduxj = (dsduxj - rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1136 >    drtduyj = (dsduyj - rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1137 >    drtduzj = (dsduzj - rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1138  
1139 + !!$    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1140 + !!$    write(*,*) 'drtdu_j = ', drtduxj, drtduyj, drtduzj
1141  
1056 –    drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1057 –    drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1058 –    drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1059 –    drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1060 –    drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1061 –    drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1062 –    drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1063 –    drtdyj = (dsdyj + rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1064 –    drtdzj = (dsdzj + rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1065 –    drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1066 –    drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1067 –    drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1068 –    
1142      rt2 = rt*rt
1143      rt3 = rt2*rt
1144      rt5 = rt2*rt3
# Line 1074 | Line 1147 | contains  
1147      rt12 = rt6*rt6
1148      rt126 = rt12 - rt6
1149  
1150 +    pot_temp = 4.0d0 * eps * rt126
1151 +
1152 +    vpair = vpair + pot_temp
1153      if (do_pot) then
1154   #ifdef IS_MPI
1155 <       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1156 <       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1155 >       pot_row(SHAPE_POT,atom1) = pot_row(SHAPE_POT,atom1) + 0.5d0*pot_temp*sw
1156 >       pot_col(SHAPE_POT,atom2) = pot_col(SHAPE_POT,atom2) + 0.5d0*pot_temp*sw
1157   #else
1158 <       pot = pot + 4.0d0*eps*rt126*sw
1158 >       pot = pot + pot_temp*sw
1159   #endif
1160      endif
1085 –    
1086 –    write(*,*) 'drtdxi = ', drtdxi, drtdyi
1087 –    write(*,*) 'depsdxi = ', depsdxi, depsdyi
1161  
1162 + !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1163 +
1164      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1165      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1166      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1093 | Line 1168 | contains  
1168      dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1169      dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1170  
1096 –    write(*,*) 'drtduzj = ', drtduzj, depsduzj
1097 –
1171      dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1172      dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1173      dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
1174      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1175      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1176      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1177 <
1177 > !!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1178      ! do the torques first since they are easy:
1179      ! remember that these are still in the body fixed axes
1180  
1181 <    write(*,*) 'dvdx = ', dvdxi, dvdyi, dvdzi
1182 <    write(*,*) 'dvdx = ', dvdxj, dvdyj, dvdzj
1183 <    write(*,*) 'dvdu = ', dvduxi, dvduyi, dvduzi
1111 <    write(*,*) 'dvdu = ', dvduxj, dvduyj, dvduzj
1181 >    txi = 0.0d0
1182 >    tyi = 0.0d0
1183 >    tzi = 0.0d0
1184  
1185 <    txi = dvduxi * sw
1186 <    tyi = dvduyi * sw
1187 <    tzi = dvduzi * sw
1185 >    txj = 0.0d0
1186 >    tyj = 0.0d0
1187 >    tzj = 0.0d0
1188  
1189 <    txj = dvduxj * sw
1190 <    tyj = dvduyj * sw
1191 <    tzj = dvduzj * sw
1189 >    txi = (dvduyi - dvduzi) * sw
1190 >    tyi = (dvduzi - dvduxi) * sw
1191 >    tzi = (dvduxi - dvduyi) * sw
1192  
1193 +    txj = (dvduyj - dvduzj) * sw
1194 +    tyj = (dvduzj - dvduxj) * sw
1195 +    tzj = (dvduxj - dvduyj) * sw
1196 +
1197 + !!$    txi = dvduxi * sw
1198 + !!$    tyi = dvduyi * sw
1199 + !!$    tzi = dvduzi * sw
1200 + !!$
1201 + !!$    txj = dvduxj * sw
1202 + !!$    tyj = dvduyj * sw
1203 + !!$    tzj = dvduzj * sw
1204 +
1205 +    write(*,*) 't1 = ', txi, tyi, tzi
1206 +    write(*,*) 't2 = ', txj, tyj, tzj
1207 +
1208      ! go back to lab frame using transpose of rotation matrix:
1209 <    
1209 >
1210   #ifdef IS_MPI
1211      t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
1212           a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
# Line 1127 | Line 1214 | contains  
1214           a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
1215      t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
1216           a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
1217 <    
1217 >
1218      t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
1219           a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
1220      t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
1221 <            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1221 >         a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1222      t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
1223           a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
1224   #else
1225      t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
1226      t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
1227      t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
1228 <    
1228 >
1229      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1230      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1231      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1232 +
1233   #endif
1234      ! Now, on to the forces:
1235 <    
1235 >
1236      ! first rotate the i terms back into the lab frame:
1149 –    
1150 –    fxi = dvdxi * sw
1151 –    fyi = dvdyi * sw
1152 –    fzi = dvdzi * sw
1237  
1238 <    fxj = dvdxj * sw
1239 <    fyj = dvdyj * sw
1240 <    fzj = dvdzj * sw
1238 >    fxi = -dvdxi * sw
1239 >    fyi = -dvdyi * sw
1240 >    fzi = -dvdzi * sw
1241  
1242 +    fxj = -dvdxj * sw
1243 +    fyj = -dvdyj * sw
1244 +    fzj = -dvdzj * sw
1245 +
1246 +
1247   #ifdef IS_MPI
1248      fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1249      fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
# Line 1167 | Line 1256 | contains  
1256      fxii = a(1,atom1)*fxi + a(4,atom1)*fyi + a(7,atom1)*fzi
1257      fyii = a(2,atom1)*fxi + a(5,atom1)*fyi + a(8,atom1)*fzi
1258      fzii = a(3,atom1)*fxi + a(6,atom1)*fyi + a(9,atom1)*fzi
1259 <    
1259 >
1260      fxjj = a(1,atom2)*fxj + a(4,atom2)*fyj + a(7,atom2)*fzj
1261      fyjj = a(2,atom2)*fxj + a(5,atom2)*fyj + a(8,atom2)*fzj
1262      fzjj = a(3,atom2)*fxj + a(6,atom2)*fyj + a(9,atom2)*fzj
# Line 1176 | Line 1265 | contains  
1265      fxij = -fxii
1266      fyij = -fyii
1267      fzij = -fzii
1268 <    
1268 >
1269      fxji = -fxjj
1270      fyji = -fyjj
1271      fzji = -fzjj
1272  
1273 <    fxradial = fxii + fxji
1274 <    fyradial = fyii + fyji
1275 <    fzradial = fzii + fzji
1276 <
1273 >    fxradial = (fxii + fxji)
1274 >    fyradial = (fyii + fyji)
1275 >    fzradial = (fzii + fzji)
1276 > !!$    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1277   #ifdef IS_MPI
1278      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1279      f_Row(2,atom1) = f_Row(2,atom1) + fyradial
1280      f_Row(3,atom1) = f_Row(3,atom1) + fzradial
1281 <    
1281 >
1282      f_Col(1,atom2) = f_Col(1,atom2) - fxradial
1283      f_Col(2,atom2) = f_Col(2,atom2) - fyradial
1284      f_Col(3,atom2) = f_Col(3,atom2) - fzradial
# Line 1197 | Line 1286 | contains  
1286      f(1,atom1) = f(1,atom1) + fxradial
1287      f(2,atom1) = f(2,atom1) + fyradial
1288      f(3,atom1) = f(3,atom1) + fzradial
1289 <    
1289 >
1290      f(1,atom2) = f(1,atom2) - fxradial
1291      f(2,atom2) = f(2,atom2) - fyradial
1292      f(3,atom2) = f(3,atom2) - fzradial
# Line 1210 | Line 1299 | contains  
1299      id1 = atom1
1300      id2 = atom2
1301   #endif
1302 <    
1302 >
1303      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1304 <      
1304 >
1305         fpair(1) = fpair(1) + fxradial
1306         fpair(2) = fpair(2) + fyradial
1307         fpair(3) = fpair(3) + fzradial
1308 <      
1308 >
1309      endif
1221 –    
1222 –    write(*,*) 'f1 = ', f(1,atom1), f(2,atom1), f(3,atom1)
1223 –    write(*,*) 'f2 = ', f(1,atom2), f(2,atom2), f(3,atom2)
1224 –    write(*,*) 't1 = ', t(1,atom1), t(2,atom1), t(3,atom1)
1225 –    write(*,*) 't2 = ', t(1,atom2), t(2,atom2), t(3,atom2)
1310  
1227 –
1311    end subroutine do_shape_pair
1312 <    
1312 >
1313    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1314  
1315      ! Purpose: Compute the associated Legendre functions
# Line 1244 | Line 1327 | contains  
1327      !
1328      ! The original Fortran77 codes can be found here:
1329      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1330 <    
1330 >
1331      real (kind=dp), intent(in) :: x
1332      integer, intent(in) :: l, m, lmax
1333      real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
# Line 1261 | Line 1344 | contains  
1344  
1345      ! start with 0,0:
1346      PLM(0,0) = 1.0D0
1347 <  
1347 >
1348      ! x = +/- 1 functions are easy:
1349      IF (abs(X).EQ.1.0D0) THEN
1350         DO I = 1, m
# Line 1316 | Line 1399 | contains  
1399  
1400  
1401    subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1402 <  
1402 >
1403      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1404      !          or Ln(x) or Hn(x), and their derivatives
1405      ! Input :  function_type --- Function code
# Line 1335 | Line 1418 | contains  
1418      !
1419      ! The original Fortran77 codes can be found here:
1420      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1421 <  
1421 >
1422      real(kind=8), intent(in) :: x
1423      integer, intent(in):: m, mmax
1424      integer, intent(in):: function_type
1425      real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1426 <  
1426 >
1427      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1428      integer :: k
1429  
# Line 1386 | Line 1469 | contains  
1469  
1470  
1471      RETURN
1472 <    
1472 >
1473    end subroutine Orthogonal_Polynomial
1391 –  
1392 – end module shapes
1474  
1475 < subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1476 <     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1396 <     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1397 <     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1398 <     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1399 <     myATID, status)
1475 >  subroutine deallocateShapes(this)
1476 >    type(Shape), pointer :: this
1477  
1478 <  use definitions
1479 <  use shapes, only: newShapeType
1480 <  
1481 <  integer :: nContactFuncs
1405 <  integer :: nRangeFuncs
1406 <  integer :: nStrengthFuncs
1407 <  integer :: status
1408 <  integer :: myATID
1409 <  
1410 <  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1411 <  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1412 <  integer, dimension(nContactFuncs) :: ContactFunctionType        
1413 <  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1414 <  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1415 <  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1416 <  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1417 <  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1418 <  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1419 <  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1420 <  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1421 <  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1422 <  
1423 <  call newShapeType(nContactFuncs, ContactFuncLValue, &
1424 <       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1425 <       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1426 <       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1427 <       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1428 <       myATID, status)
1478 >    if (associated( this%ContactFuncLValue)) then
1479 >       deallocate(this%ContactFuncLValue)
1480 >       this%ContactFuncLValue => null()
1481 >    end if
1482  
1483 <  return
1484 < end subroutine makeShape
1483 >    if (associated( this%ContactFuncMValue)) then
1484 >       deallocate( this%ContactFuncMValue)
1485 >       this%ContactFuncMValue => null()
1486 >    end if
1487 >    if (associated( this%ContactFunctionType)) then
1488 >       deallocate(this%ContactFunctionType)
1489 >       this%ContactFunctionType => null()
1490 >    end if
1491  
1492 < subroutine completeShapeFF(status)
1492 >    if (associated( this%ContactFuncCoefficient)) then
1493 >       deallocate(this%ContactFuncCoefficient)
1494 >       this%ContactFuncCoefficient => null()
1495 >    end if
1496  
1497 <  use shapes, only: complete_Shape_FF
1497 >    if (associated( this%RangeFuncLValue)) then
1498 >       deallocate(this%RangeFuncLValue)
1499 >       this%RangeFuncLValue => null()
1500 >    end if
1501 >    if (associated( this%RangeFuncMValue)) then
1502 >       deallocate( this%RangeFuncMValue)
1503 >       this%RangeFuncMValue => null()
1504 >    end if
1505  
1506 <  integer, intent(out)  :: status
1507 <  integer :: myStatus
1506 >    if (associated( this%RangeFunctionType)) then
1507 >       deallocate( this%RangeFunctionType)
1508 >       this%RangeFunctionType => null()
1509 >    end if
1510 >    if (associated( this%RangeFuncCoefficient)) then
1511 >       deallocate(this%RangeFuncCoefficient)
1512 >       this%RangeFuncCoefficient => null()
1513 >    end if
1514  
1515 <  myStatus = 0
1515 >    if (associated( this%StrengthFuncLValue)) then
1516 >       deallocate(this%StrengthFuncLValue)
1517 >       this%StrengthFuncLValue => null()
1518 >    end if
1519  
1520 <  call complete_Shape_FF(myStatus)
1520 >    if (associated( this%StrengthFuncMValue )) then
1521 >       deallocate(this%StrengthFuncMValue)
1522 >       this%StrengthFuncMValue => null()
1523 >    end if
1524  
1525 <  status = myStatus
1525 >    if(associated( this%StrengthFunctionType)) then
1526 >       deallocate(this%StrengthFunctionType)
1527 >       this%StrengthFunctionType => null()
1528 >    end if
1529 >    if (associated( this%StrengthFuncCoefficient )) then
1530 >       deallocate(this%StrengthFuncCoefficient)
1531 >       this%StrengthFuncCoefficient => null()
1532 >    end if
1533 >  end subroutine deallocateShapes
1534  
1535 <  return
1536 < end subroutine completeShapeFF
1535 >  subroutine destroyShapeTypes
1536 >    integer :: i
1537 >    type(Shape), pointer :: thisShape
1538  
1539 +    ! First walk through and kill the shape
1540 +    do i = 1,ShapeMap%n_shapes
1541 +       thisShape => ShapeMap%Shapes(i)
1542 +       call deallocateShapes(thisShape)
1543 +    end do
1544 +
1545 +    ! set shape map to starting values
1546 +    ShapeMap%n_shapes = 0
1547 +    ShapeMap%currentShape = 0
1548 +
1549 +    if (associated(ShapeMap%Shapes)) then
1550 +       deallocate(ShapeMap%Shapes)
1551 +       ShapeMap%Shapes => null()
1552 +    end if
1553 +
1554 +    if (associated(ShapeMap%atidToShape)) then
1555 +       deallocate(ShapeMap%atidToShape)
1556 +       ShapeMap%atidToShape => null()
1557 +    end if
1558 +
1559 +
1560 +  end subroutine destroyShapeTypes
1561 +
1562 +
1563 + end module shapes

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines