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 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 55 | 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 68 | Line 69 | module shapes
69    public :: do_shape_pair
70    public :: newShapeType
71    public :: complete_Shape_FF
72 +  public :: destroyShapeTypes
73 +  public :: getShapeCut
74  
72 –
75    type, private :: Shape
76       integer :: atid
77       integer :: nContactFuncs
# Line 93 | 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 144 | 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 167 | 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 189 | 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 227 | 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 331 | 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 345 | 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
362 –          
366            ShapeMap%currentShape = ShapeMap%currentShape + 1
367            current = ShapeMap%currentShape
368  
369 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
370 <          ShapeMap%atidToShape(thisIP) = current
368 <          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 453 | 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 480 | 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 489 | Line 509 | contains  
509         call handleError("calc_shape", "NO SHAPEMAP!!!!")
510         return      
511      endif
512 <    
512 >
513      !! We assume that the rotation matrices have already been calculated
514      !! and placed in the A array.
495 –
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.0_dp
522 +    drduyi = 0.0_dp
523 +    drduzi = 0.0_dp
524  
525      drdxj = d(1) / rij
526      drdyj = d(2) / rij
527      drdzj = d(3) / rij
528 <    
528 >    drduxj = 0.0_dp
529 >    drduyj = 0.0_dp
530 >    drduzj = 0.0_dp
531 >
532      ! find the atom type id (atid) for each atom:
533   #ifdef IS_MPI
534      atid1 = atid_Row(atom1)
# Line 516 | 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  
549         sigma_i = ShapeMap%Shapes(st1)%sigma
550         s_i = ShapeMap%Shapes(st1)%sigma
551         eps_i = ShapeMap%Shapes(st1)%epsilon
552 <       dsigmaidx = 0.0d0
553 <       dsigmaidy = 0.0d0
554 <       dsigmaidz = 0.0d0
555 <       dsigmaidux = 0.0d0
556 <       dsigmaiduy = 0.0d0
557 <       dsigmaiduz = 0.0d0
558 <       dsidx = 0.0d0
559 <       dsidy = 0.0d0
560 <       dsidz = 0.0d0
561 <       dsidux = 0.0d0
562 <       dsiduy = 0.0d0
563 <       dsiduz = 0.0d0
564 <       depsidx = 0.0d0
565 <       depsidy = 0.0d0
566 <       depsidz = 0.0d0
567 <       depsidux = 0.0d0
568 <       depsiduy = 0.0d0
569 <       depsiduz = 0.0d0
552 >       dsigmaidx = 0.0_dp
553 >       dsigmaidy = 0.0_dp
554 >       dsigmaidz = 0.0_dp
555 >       dsigmaidux = 0.0_dp
556 >       dsigmaiduy = 0.0_dp
557 >       dsigmaiduz = 0.0_dp
558 >       dsidx = 0.0_dp
559 >       dsidy = 0.0_dp
560 >       dsidz = 0.0_dp
561 >       dsidux = 0.0_dp
562 >       dsiduy = 0.0_dp
563 >       dsiduz = 0.0_dp
564 >       depsidx = 0.0_dp
565 >       depsidy = 0.0_dp
566 >       depsidz = 0.0_dp
567 >       depsidux = 0.0_dp
568 >       depsiduy = 0.0_dp
569 >       depsiduz = 0.0_dp
570      else
571  
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)
560 –      
561 – #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 570 | Line 599 | contains  
599  
600         dctidx = - zi * xi / r3
601         dctidy = - zi * yi / r3
602 <       dctidz = 1.0d0 / rij - zi2 / r3
603 <       dctidux = - (zi * xi2) / r3
604 <       dctiduy = - (zi * yi2) / r3
605 <       dctiduz = zi / rij - (zi2 * zi) / r3
602 >       dctidz = 1.0_dp / rij - zi2 / r3
603 >       dctidux = yi / rij ! - (zi * xi2) / r3
604 >       dctiduy = -xi / rij !- (zi * yi2) / r3
605 >       dctiduz = 0.0_dp !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:
609  
610         sti2 = 1.0_dp - cti*cti
611 <       if (dabs(sti2) .lt. 1.0d-12) then
612 <          proji = sqrt(rij * 1.0d-12)
613 <          dcpidx = 1.0d0 / proji
614 <          dcpidy = 0.0d0
611 >       if (abs(sti2) .lt. 1.0e-12_dp) then
612 >          proji = sqrt(rij * 1.0e-12_dp)
613 >          dcpidx = 1.0_dp / proji
614 >          dcpidy = 0.0_dp
615            dcpidux = xi / proji
616 <          dcpiduy = 0.0d0
617 <          dspidx = 0.0d0
618 <          dspidy = 1.0d0 / proji
619 <          dspidux = 0.0d0
616 >          dcpiduy = 0.0_dp
617 >          dspidx = 0.0_dp
618 >          dspidy = 1.0_dp / proji
619 >          dspidux = 0.0_dp
620            dspiduy = yi / proji
621         else
622            proji = sqrt(xi2 + yi2)
623            proji3 = proji*proji*proji
624 <          dcpidx = 1.0d0 / proji - xi2 / proji3
624 >          dcpidx = 1.0_dp / proji - xi2 / proji3
625            dcpidy = - xi * yi / 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
629 >          dspidy = 1.0_dp / proji - yi2 / 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 = 0.0d0
637 <      
635 >       dcpidz = 0.0_dp
636 >       dcpiduz = 0.0_dp
637 >
638         spi = yi / proji
639 <       dspidz = 0.0d0
640 <       dspiduz = 0.0d0
639 >       dspidz = 0.0_dp
640 >       dspiduz = 0.0_dp
641  
642         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
643              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 618 | 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)
621 –      
622 –       sigma_i = 0.0d0
623 –       s_i = 0.0d0
624 –       eps_i = 0.0d0
625 –       dsigmaidx = 0.0d0
626 –       dsigmaidy = 0.0d0
627 –       dsigmaidz = 0.0d0
628 –       dsigmaidux = 0.0d0
629 –       dsigmaiduy = 0.0d0
630 –       dsigmaiduz = 0.0d0
631 –       dsidx = 0.0d0
632 –       dsidy = 0.0d0
633 –       dsidz = 0.0d0
634 –       dsidux = 0.0d0
635 –       dsiduy = 0.0d0
636 –       dsiduz = 0.0d0
637 –       depsidx = 0.0d0
638 –       depsidy = 0.0d0
639 –       depsidz = 0.0d0
640 –       depsidux = 0.0d0
641 –       depsiduy = 0.0d0
642 –       depsiduz = 0.0d0
650  
651 +       sigma_i = 0.0_dp
652 +       s_i = 0.0_dp
653 +       eps_i = 0.0_dp
654 +       dsigmaidx = 0.0_dp
655 +       dsigmaidy = 0.0_dp
656 +       dsigmaidz = 0.0_dp
657 +       dsigmaidux = 0.0_dp
658 +       dsigmaiduy = 0.0_dp
659 +       dsigmaiduz = 0.0_dp
660 +       dsidx = 0.0_dp
661 +       dsidy = 0.0_dp
662 +       dsidz = 0.0_dp
663 +       dsidux = 0.0_dp
664 +       dsiduy = 0.0_dp
665 +       dsiduz = 0.0_dp
666 +       depsidx = 0.0_dp
667 +       depsidy = 0.0_dp
668 +       depsidz = 0.0_dp
669 +       depsidux = 0.0_dp
670 +       depsiduy = 0.0_dp
671 +       depsiduz = 0.0_dp
672 +
673         do lm = 1, ShapeMap%Shapes(st1)%nContactFuncs
674            l = ShapeMap%Shapes(st1)%ContactFuncLValue(lm)
675            m = ShapeMap%Shapes(st1)%ContactFuncMValue(lm)
# Line 652 | 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 666 | 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
676 –          
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 688 | 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)
722 <          
722 >
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 708 | Line 739 | contains  
739            endif
740  
741            s_i = s_i + plm_i(m,l)*Phunc
742 <          
742 >
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 724 | 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 750 | 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 768 | Line 799 | contains  
799         end do
800  
801      endif
771 –      
772 –       ! 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
808         eps_j = ShapeMap%Shapes(st2)%epsilon
809 <       dsigmajdx = 0.0d0
810 <       dsigmajdy = 0.0d0
811 <       dsigmajdz = 0.0d0
812 <       dsigmajdux = 0.0d0
813 <       dsigmajduy = 0.0d0
814 <       dsigmajduz = 0.0d0
815 <       dsjdx = 0.0d0
816 <       dsjdy = 0.0d0
817 <       dsjdz = 0.0d0
818 <       dsjdux = 0.0d0
819 <       dsjduy = 0.0d0
820 <       dsjduz = 0.0d0
821 <       depsjdx = 0.0d0
822 <       depsjdy = 0.0d0
823 <       depsjdz = 0.0d0
824 <       depsjdux = 0.0d0
825 <       depsjduy = 0.0d0
826 <       depsjduz = 0.0d0
809 >       dsigmajdx = 0.0_dp
810 >       dsigmajdy = 0.0_dp
811 >       dsigmajdz = 0.0_dp
812 >       dsigmajdux = 0.0_dp
813 >       dsigmajduy = 0.0_dp
814 >       dsigmajduz = 0.0_dp
815 >       dsjdx = 0.0_dp
816 >       dsjdy = 0.0_dp
817 >       dsjdz = 0.0_dp
818 >       dsjdux = 0.0_dp
819 >       dsjduy = 0.0_dp
820 >       dsjduz = 0.0_dp
821 >       depsjdx = 0.0_dp
822 >       depsjdy = 0.0_dp
823 >       depsjdz = 0.0_dp
824 >       depsjdux = 0.0_dp
825 >       depsjduy = 0.0_dp
826 >       depsjduz = 0.0_dp
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 807 | 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 = - (zi * xj2) / r3
862 <       dctjduy = - (zj * yj2) / r3
863 <       dctjduz = zj / rij - (zj2 * zj) / r3
864 <      
860 >       dctjdz = 1.0_dp / rij - zj2 / r3
861 >       dctjdux = yj / rij !- (zi * xj2) / r3
862 >       dctjduy = -xj / rij !- (zj * yj2) / r3
863 >       dctjduz = 0.0_dp !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  
868         stj2 = 1.0_dp - ctj*ctj
869 <       if (dabs(stj2) .lt. 1.0d-12) then
870 <          projj = sqrt(rij * 1.0d-12)
871 <          dcpjdx = 1.0d0 / projj
872 <          dcpjdy = 0.0d0
869 >       if (abs(stj2) .lt. 1.0e-12_dp) then
870 >          projj = sqrt(rij * 1.0e-12_dp)
871 >          dcpjdx = 1.0_dp / projj
872 >          dcpjdy = 0.0_dp
873            dcpjdux = xj / projj
874 <          dcpjduy = 0.0d0
875 <          dspjdx = 0.0d0
876 <          dspjdy = 1.0d0 / projj
877 <          dspjdux = 0.0d0
874 >          dcpjduy = 0.0_dp
875 >          dspjdx = 0.0_dp
876 >          dspjdy = 1.0_dp / projj
877 >          dspjdux = 0.0_dp
878            dspjduy = yj / projj
879         else
880            projj = sqrt(xj2 + yj2)
881            projj3 = projj*projj*projj
882 <          dcpjdx = 1.0d0 / projj - xj2 / projj3
882 >          dcpjdx = 1.0_dp / projj - xj2 / projj3
883            dcpjdy = - xj * yj / 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
887 >          dspjdy = 1.0_dp / projj - yj2 / 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 = 0.0d0
895 <      
893 >       dcpjdz = 0.0_dp
894 >       dcpjduz = 0.0_dp
895 >
896         spj = yj / projj
897 <       dspjdz = 0.0d0
898 <       dspjduz = 0.0d0
897 >       dspjdz = 0.0_dp
898 >       dspjduz = 0.0_dp
899  
900  
901 <       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
902 <       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
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)
877 –      
878 –       sigma_j = 0.0d0
879 –       s_j = 0.0d0
880 –       eps_j = 0.0d0
881 –       dsigmajdx = 0.0d0
882 –       dsigmajdy = 0.0d0
883 –       dsigmajdz = 0.0d0
884 –       dsigmajdux = 0.0d0
885 –       dsigmajduy = 0.0d0
886 –       dsigmajduz = 0.0d0
887 –       dsjdx = 0.0d0
888 –       dsjdy = 0.0d0
889 –       dsjdz = 0.0d0
890 –       dsjdux = 0.0d0
891 –       dsjduy = 0.0d0
892 –       dsjduz = 0.0d0
893 –       depsjdx = 0.0d0
894 –       depsjdy = 0.0d0
895 –       depsjdz = 0.0d0
896 –       depsjdux = 0.0d0
897 –       depsjduy = 0.0d0
898 –       depsjduz = 0.0d0
911  
912 +       sigma_j = 0.0_dp
913 +       s_j = 0.0_dp
914 +       eps_j = 0.0_dp
915 +       dsigmajdx = 0.0_dp
916 +       dsigmajdy = 0.0_dp
917 +       dsigmajdz = 0.0_dp
918 +       dsigmajdux = 0.0_dp
919 +       dsigmajduy = 0.0_dp
920 +       dsigmajduz = 0.0_dp
921 +       dsjdx = 0.0_dp
922 +       dsjdy = 0.0_dp
923 +       dsjdz = 0.0_dp
924 +       dsjdux = 0.0_dp
925 +       dsjduy = 0.0_dp
926 +       dsjduz = 0.0_dp
927 +       depsjdx = 0.0_dp
928 +       depsjdy = 0.0_dp
929 +       depsjdz = 0.0_dp
930 +       depsjdux = 0.0_dp
931 +       depsjduy = 0.0_dp
932 +       depsjduz = 0.0_dp
933 +
934         do lm = 1, ShapeMap%Shapes(st2)%nContactFuncs
935            l = ShapeMap%Shapes(st2)%ContactFuncLValue(lm)
936            m = ShapeMap%Shapes(st2)%ContactFuncMValue(lm)
# Line 908 | 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 920 | 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 950 | 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 964 | 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 1005 | 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
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 1030 | 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 1062 | Line 1096 | contains  
1096      dsduzj = 0.5*dsjduz
1097  
1098      eps = sqrt(eps_i * eps_j)
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.0_dp * eps)
1103 +    depsdyi = eps_j * depsidy / (2.0_dp * eps)
1104 +    depsdzi = eps_j * depsidz / (2.0_dp * eps)
1105 +    depsduxi = eps_j * depsidux / (2.0_dp * eps)
1106 +    depsduyi = eps_j * depsiduy / (2.0_dp * eps)
1107 +    depsduzi = eps_j * depsiduz / (2.0_dp * eps)
1108  
1109 <    depsdxi = eps_j * depsidx / (2.0d0 * eps)
1110 <    depsdyi = eps_j * depsidy / (2.0d0 * eps)
1111 <    depsdzi = eps_j * depsidz / (2.0d0 * eps)
1112 <    depsduxi = eps_j * depsidux / (2.0d0 * eps)
1113 <    depsduyi = eps_j * depsiduy / (2.0d0 * eps)
1114 <    depsduzi = eps_j * depsiduz / (2.0d0 * eps)
1109 >    depsdxj = eps_i * depsjdx / (2.0_dp * eps)
1110 >    depsdyj = eps_i * depsjdy / (2.0_dp * eps)
1111 >    depsdzj = eps_i * depsjdz / (2.0_dp * eps)
1112 >    depsduxj = eps_i * depsjdux / (2.0_dp * eps)
1113 >    depsduyj = eps_i * depsjduy / (2.0_dp * eps)
1114 >    depsduzj = eps_i * depsjduz / (2.0_dp * eps)
1115  
1073 –    depsdxj = eps_i * depsjdx / (2.0d0 * eps)
1074 –    depsdyj = eps_i * depsjdy / (2.0d0 * eps)
1075 –    depsdzj = eps_i * depsjdz / (2.0d0 * eps)
1076 –    depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1077 –    depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1078 –    depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1079 –    
1116   !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1117 +
1118   !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1082 – !!$
1083 – !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1119   !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1120   !!$
1121   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
# Line 1088 | Line 1123 | contains  
1123      rtdenom = rij-sigma+s
1124      rt = s / rtdenom
1125  
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 <    
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 >
1142      rt2 = rt*rt
1143      rt3 = rt2*rt
1144      rt5 = rt2*rt3
# Line 1109 | Line 1147 | contains  
1147      rt12 = rt6*rt6
1148      rt126 = rt12 - rt6
1149  
1150 <    pot_temp = 4.0d0 * eps * rt126
1150 >    pot_temp = 4.0_dp * eps * rt126
1151  
1152      vpair = vpair + pot_temp
1153      if (do_pot) then
1154   #ifdef IS_MPI
1155 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1156 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1155 >       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5_dp*pot_temp*sw
1156 >       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5_dp*pot_temp*sw
1157   #else
1158         pot = pot + pot_temp*sw
1159   #endif
1160      endif
1161  
1162   !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1125 –    
1126 –    dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1127 –    dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1128 –    dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
1129 –    dvduxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxi + 4.0d0*depsduxi*rt126
1130 –    dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1131 –    dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1163  
1164 <    dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1165 <    dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1166 <    dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
1167 <    dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1168 <    dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1169 <    dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1164 >    dvdxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxi + 4.0_dp*depsdxi*rt126
1165 >    dvdyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyi + 4.0_dp*depsdyi*rt126
1166 >    dvdzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzi + 4.0_dp*depsdzi*rt126
1167 >    dvduxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxi + 4.0_dp*depsduxi*rt126
1168 >    dvduyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyi + 4.0_dp*depsduyi*rt126
1169 >    dvduzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzi + 4.0_dp*depsduzi*rt126
1170  
1171 +    dvdxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxj + 4.0_dp*depsdxj*rt126
1172 +    dvdyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyj + 4.0_dp*depsdyj*rt126
1173 +    dvdzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzj + 4.0_dp*depsdzj*rt126
1174 +    dvduxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxj + 4.0_dp*depsduxj*rt126
1175 +    dvduyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyj + 4.0_dp*depsduyj*rt126
1176 +    dvduzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzj + 4.0_dp*depsduzj*rt126
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 +    txi = 0.0_dp
1182 +    tyi = 0.0_dp
1183 +    tzi = 0.0_dp
1184  
1185 < !!$    write(*,*) 'sw = ', sw
1186 < !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1187 < !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1147 < !!$
1148 <    txi =  (dvduzi - dvduyi) * sw
1149 <    tyi =  (dvduxi - dvduzi) * sw
1150 <    tzi =  (dvduyi - dvduxi) * sw
1185 >    txj = 0.0_dp
1186 >    tyj = 0.0_dp
1187 >    tzj = 0.0_dp
1188  
1189 <    txj = (dvduzj - dvduyj) * sw
1190 <    tyj = (dvduxj - dvduzj) * sw
1191 <    tzj = (dvduyj - dvduxj) * sw
1189 >    txi = (dvduyi - dvduzi) * sw
1190 >    tyi = (dvduzi - dvduxi) * sw
1191 >    tzi = (dvduxi - dvduyi) * sw
1192  
1193 < !!$    txi = -dvduxi * sw
1194 < !!$    tyi = -dvduyi * sw
1195 < !!$    tzi = -dvduzi * sw
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 <
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 1173 | 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:
1195 –    
1196 –    fxi = dvdxi * sw
1197 –    fyi = dvdyi * sw
1198 –    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 1213 | 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 1222 | 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 = 0.5_dp * (fxii + fxji)
1274 <    fyradial = 0.5_dp * (fyii + fyji)
1275 <    fzradial = 0.5_dp * (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 1243 | 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 1256 | 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
1310  
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 1284 | 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 1300 | Line 1343 | contains  
1343      end DO
1344  
1345      ! start with 0,0:
1346 <    PLM(0,0) = 1.0D0
1347 <  
1346 >    PLM(0,0) = 1.0_DP
1347 >
1348      ! x = +/- 1 functions are easy:
1349 <    IF (abs(X).EQ.1.0D0) THEN
1349 >    IF (abs(X).EQ.1.0_DP) THEN
1350         DO I = 1, m
1351            PLM(0, I) = X**I
1352 <          DLM(0, I) = 0.5D0*I*(I+1.0D0)*X**(I+1)
1352 >          DLM(0, I) = 0.5_DP*I*(I+1.0_DP)*X**(I+1)
1353         end DO
1354         DO J = 1, m
1355            DO I = 1, l
1356               IF (I.EQ.1) THEN
1357                  DLM(I, J) = 1.0D+300
1358               ELSE IF (I.EQ.2) THEN
1359 <                DLM(I, J) = -0.25D0*(J+2)*(J+1)*J*(J-1)*X**(J+1)
1359 >                DLM(I, J) = -0.25_DP*(J+2)*(J+1)*J*(J-1)*X**(J+1)
1360               ENDIF
1361            end DO
1362         end DO
# Line 1321 | Line 1364 | contains  
1364      ENDIF
1365  
1366      LS = 1
1367 <    IF (abs(X).GT.1.0D0) LS = -1
1368 <    XQ = sqrt(LS*(1.0D0-X*X))
1369 <    XS = LS*(1.0D0-X*X)
1367 >    IF (abs(X).GT.1.0_DP) LS = -1
1368 >    XQ = sqrt(LS*(1.0_DP-X*X))
1369 >    XS = LS*(1.0_DP-X*X)
1370  
1371      DO I = 1, l
1372 <       PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1372 >       PLM(I, I) = -LS*(2.0_DP*I-1.0_DP)*XQ*PLM(I-1, I-1)
1373      enddo
1374  
1375      DO I = 0, l
1376 <       PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1376 >       PLM(I, I+1)=(2.0_DP*I+1.0_DP)*X*PLM(I, I)
1377      enddo
1378  
1379      DO I = 0, l
1380         DO J = I+2, m
1381 <          PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1382 <               (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1381 >          PLM(I, J)=((2.0_DP*J-1.0_DP)*X*PLM(I,J-1) - &
1382 >               (I+J-1.0_DP)*PLM(I,J-2))/(J-I)
1383         end DO
1384      end DO
1385  
1386 <    DLM(0, 0)=0.0D0
1386 >    DLM(0, 0)=0.0_DP
1387      DO J = 1, m
1388         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1389      end DO
1390  
1391      DO I = 1, l
1392         DO J = I, m
1393 <          DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1393 >          DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0_DP)/XQ*PLM(I-1, J)
1394         end DO
1395      end DO
1396  
# Line 1356 | 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 1375 | 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  
1430 <    A = 2.0D0
1431 <    B = 0.0D0
1432 <    C = 1.0D0
1433 <    Y0 = 1.0D0
1434 <    Y1 = 2.0D0*X
1435 <    DY0 = 0.0D0
1436 <    DY1 = 2.0D0
1437 <    PL(0) = 1.0D0
1438 <    PL(1) = 2.0D0*X
1439 <    DPL(0) = 0.0D0
1440 <    DPL(1) = 2.0D0
1430 >    A = 2.0_DP
1431 >    B = 0.0_DP
1432 >    C = 1.0_DP
1433 >    Y0 = 1.0_DP
1434 >    Y1 = 2.0_DP*X
1435 >    DY0 = 0.0_DP
1436 >    DY1 = 2.0_DP
1437 >    PL(0) = 1.0_DP
1438 >    PL(1) = 2.0_DP*X
1439 >    DPL(0) = 0.0_DP
1440 >    DPL(1) = 2.0_DP
1441      IF (function_type.EQ.CHEBYSHEV_TN) THEN
1442         Y1 = X
1443 <       DY1 = 1.0D0
1443 >       DY1 = 1.0_DP
1444         PL(1) = X
1445 <       DPL(1) = 1.0D0
1445 >       DPL(1) = 1.0_DP
1446      ELSE IF (function_type.EQ.LAGUERRE) THEN
1447 <       Y1 = 1.0D0-X
1448 <       DY1 = -1.0D0
1449 <       PL(1) = 1.0D0-X
1450 <       DPL(1) = -1.0D0
1447 >       Y1 = 1.0_DP-X
1448 >       DY1 = -1.0_DP
1449 >       PL(1) = 1.0_DP-X
1450 >       DPL(1) = -1.0_DP
1451      ENDIF
1452      DO K = 2, m
1453         IF (function_type.EQ.LAGUERRE) THEN
1454 <          A = -1.0D0/K
1455 <          B = 2.0D0+A
1456 <          C = 1.0D0+A
1454 >          A = -1.0_DP/K
1455 >          B = 2.0_DP+A
1456 >          C = 1.0_DP+A
1457         ELSE IF (function_type.EQ.HERMITE) THEN
1458 <          C = 2.0D0*(K-1.0D0)
1458 >          C = 2.0_DP*(K-1.0_DP)
1459         ENDIF
1460         YN = (A*X+B)*Y1-C*Y0
1461         DYN = A*Y1+(A*X+B)*DY1-C*DY0
# Line 1426 | Line 1469 | contains  
1469  
1470  
1471      RETURN
1472 <    
1472 >
1473    end subroutine Orthogonal_Polynomial
1474 <  
1474 >
1475 >  subroutine deallocateShapes(this)
1476 >    type(Shape), pointer :: this
1477 >
1478 >    if (associated( this%ContactFuncLValue)) then
1479 >       deallocate(this%ContactFuncLValue)
1480 >       this%ContactFuncLValue => null()
1481 >    end if
1482 >
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 >    if (associated( this%ContactFuncCoefficient)) then
1493 >       deallocate(this%ContactFuncCoefficient)
1494 >       this%ContactFuncCoefficient => null()
1495 >    end if
1496 >
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 >    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 >    if (associated( this%StrengthFuncLValue)) then
1516 >       deallocate(this%StrengthFuncLValue)
1517 >       this%StrengthFuncLValue => null()
1518 >    end if
1519 >
1520 >    if (associated( this%StrengthFuncMValue )) then
1521 >       deallocate(this%StrengthFuncMValue)
1522 >       this%StrengthFuncMValue => null()
1523 >    end if
1524 >
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 >  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