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 1698 by chrisfen, Tue Nov 2 15:28:25 2004 UTC vs.
Revision 2214 by chrisfen, Wed Apr 27 20:14:03 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 >
59    INTEGER, PARAMETER:: CHEBYSHEV_TN = 1
60    INTEGER, PARAMETER:: CHEBYSHEV_UN = 2
61    INTEGER, PARAMETER:: LAGUERRE     = 3
# Line 26 | Line 68 | module shapes
68    public :: do_shape_pair
69    public :: newShapeType
70    public :: complete_Shape_FF
71 +  public :: destroyShapeTypes
72  
30 –
73    type, private :: Shape
74       integer :: atid
75       integer :: nContactFuncs
# Line 51 | Line 93 | module shapes
93       real ( kind = dp )  :: epsilon
94       real ( kind = dp )  :: sigma
95    end type Shape
96 <  
96 >
97    type, private :: ShapeList
98       integer :: n_shapes = 0
99       integer :: currentShape = 0
100 <     type (Shape), pointer :: Shapes(:)      => null()
101 <     integer, pointer      :: atidToShape(:) => null()
100 >     type(Shape), pointer :: Shapes(:)      => null()
101 >     integer, pointer     :: atidToShape(:) => null()
102    end type ShapeList
103    
104    type(ShapeList), save :: ShapeMap
105 <
105 >  
106    integer :: lmax
107 <
107 >  
108   contains  
109 <
109 >  
110    subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
111         ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115 <       myATID, status)
116 <
115 >       c_ident, status)
116 >    
117      integer :: nContactFuncs
118      integer :: nRangeFuncs
119      integer :: nStrengthFuncs
120      integer :: shape_ident
121      integer :: status
122 +    integer :: c_ident
123      integer :: myATID
124      integer :: bigL
125      integer :: bigM
# Line 102 | Line 145 | contains  
145  
146         call getMatchingElementList(atypes, "is_Shape", .true., &
147              nShapeTypes, MatchList)
148 <      
148 >
149         call getMatchingElementList(atypes, "is_LennardJones", .true., &
150              nLJTypes, MatchList)
151 <      
151 >
152         ShapeMap%n_shapes = nShapeTypes + nLJTypes
153 <      
153 >
154         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
155 <      
155 >
156         ntypes = getSize(atypes)
157 <      
158 <       allocate(ShapeMap%atidToShape(0:ntypes))
157 >
158 >       allocate(ShapeMap%atidToShape(ntypes))
159      end if
160 <    
160 >
161      ShapeMap%currentShape = ShapeMap%currentShape + 1
162      current = ShapeMap%currentShape
163  
# Line 125 | Line 168 | contains  
168         return
169      endif
170  
171 <    call getElementProperty(atypes, myATID, 'c_ident', me)
171 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
172  
173 <    ShapeMap%atidToShape(me)                         = current
174 <    ShapeMap%Shapes(current)%atid                    = me
173 >    ShapeMap%atidToShape(myATID)                     = current
174 >    ShapeMap%Shapes(current)%atid                    = myATID
175      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
176      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
177      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 147 | Line 190 | contains  
190  
191      bigL = -1
192      bigM = -1
193 <    
193 >
194      do j = 1, ShapeMap%Shapes(current)%nContactFuncs
195         if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
196            bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
# Line 185 | Line 228 | contains  
228      type(Shape), intent(inout) :: myShape
229      integer, intent(out) :: stat
230      integer :: alloc_stat
231 <
231 >
232      stat = 0
233      if (associated(myShape%contactFuncLValue)) then
234         deallocate(myShape%contactFuncLValue)
# Line 289 | Line 332 | contains  
332      return
333  
334    end subroutine allocateShape
335 <    
335 >
336    subroutine complete_Shape_FF(status)
337      integer :: status
338      integer :: i, j, l, m, lm, function_type
339      real(kind=dp) :: thisDP, sigma
340 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
340 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
341      logical :: thisProperty
342  
343      status = 0
# Line 303 | Line 346 | contains  
346         status = -1
347         return
348      end if
349 <    
349 >
350      nAtypes = getSize(atypes)
351  
352      if (nAtypes == 0) then
353         status = -1
354         return
355 <    end if
355 >    end if      
356  
357      ! atypes comes from c side
358 <    do i = 0, nAtypes
359 <      
360 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
361 <      
358 >    do i = 1, nAtypes
359 >    
360 >       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
361 >       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
362 >        
363         if (thisProperty) then
320 –          
364            ShapeMap%currentShape = ShapeMap%currentShape + 1
365            current = ShapeMap%currentShape
366  
367 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
368 <          ShapeMap%atidToShape(thisIP) = current
326 <          ShapeMap%Shapes(current)%atid = thisIP
367 >          ShapeMap%atidToShape(myATID) = current
368 >          ShapeMap%Shapes(current)%atid = myATID
369  
370            ShapeMap%Shapes(current)%isLJ = .true.
371  
372 <          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
373 <          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
374 <          
372 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
373 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
374 >
375         endif
376 <      
376 >
377      end do
378  
379      haveShapeMap = .true.
380 <    
380 >
381 > !    do i = 1, ShapeMap%n_shapes
382 > !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
383 > !    end do
384 >
385    end subroutine complete_Shape_FF
386 <    
386 >
387    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
388         pot, A, f, t, do_pot)
389 <    
389 >
390      INTEGER, PARAMETER:: LMAX         = 64
391      INTEGER, PARAMETER:: MMAX         = 64
392  
# Line 348 | Line 394 | contains  
394      real (kind=dp), intent(inout) :: rij, r2
395      real (kind=dp), dimension(3), intent(in) :: d
396      real (kind=dp), dimension(3), intent(inout) :: fpair
397 <    real (kind=dp) :: pot, vpair, sw
397 >    real (kind=dp) :: pot, vpair, sw, dswdr
398      real (kind=dp), dimension(9,nLocal) :: A
399      real (kind=dp), dimension(3,nLocal) :: f
400      real (kind=dp), dimension(3,nLocal) :: t
# Line 359 | Line 405 | contains  
405      integer :: l, m, lm, id1, id2, localError, function_type
406      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
407      real (kind=dp) :: coeff
408 +    real (kind=dp) :: pot_temp
409  
410      real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
411      real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
# Line 377 | Line 424 | contains  
424  
425      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
426  
427 +    real (kind=dp) :: sti2, stj2
428 +
429      real (kind=dp) :: proji, proji3, projj, projj3
430      real (kind=dp) :: cti, ctj, cpi, cpj, spi, spj
431      real (kind=dp) :: Phunc, sigma, s, eps, rtdenom, rt
# Line 408 | Line 457 | contains  
457      real (kind=dp) :: dsduxi, dsduyi, dsduzi
458      real (kind=dp) :: dsdxj, dsdyj, dsdzj
459      real (kind=dp) :: dsduxj, dsduyj, dsduzj
460 <    
460 >
461      real (kind=dp) :: depsdxi, depsdyi, depsdzi
462      real (kind=dp) :: depsduxi, depsduyi, depsduzi
463      real (kind=dp) :: depsdxj, depsdyj, depsdzj
# Line 435 | Line 484 | contains  
484      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
485      real (kind=dp) :: fxradial, fyradial, fzradial
486  
487 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
488 +
489      real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
490      real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
491      real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
# Line 444 | Line 495 | contains  
495         call handleError("calc_shape", "NO SHAPEMAP!!!!")
496         return      
497      endif
498 <    
498 >
499      !! We assume that the rotation matrices have already been calculated
500      !! and placed in the A array.
450 –
501      r3 = r2*rij
502      r5 = r3*r2
503 <    
503 >
504      drdxi = -d(1) / rij
505      drdyi = -d(2) / rij
506      drdzi = -d(3) / rij
507 +    drduxi = 0.0d0
508 +    drduyi = 0.0d0
509 +    drduzi = 0.0d0
510  
511      drdxj = d(1) / rij
512      drdyj = d(2) / rij
513      drdzj = d(3) / rij
514 <    
514 >    drduxj = 0.0d0
515 >    drduyj = 0.0d0
516 >    drduzj = 0.0d0
517 >
518      ! find the atom type id (atid) for each atom:
519   #ifdef IS_MPI
520      atid1 = atid_Row(atom1)
# Line 471 | Line 527 | contains  
527      ! use the atid to find the shape type (st) for each atom:
528      st1 = ShapeMap%atidToShape(atid1)
529      st2 = ShapeMap%atidToShape(atid2)
530 +    
531 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
532  
533      if (ShapeMap%Shapes(st1)%isLJ) then
534  
# Line 500 | Line 558 | contains  
558   #ifdef IS_MPI
559         ! rotate the inter-particle separation into the two different
560         ! body-fixed coordinate systems:
561 <      
561 >
562         xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
563         yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
564         zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
565 <      
565 >
566   #else
567         ! rotate the inter-particle separation into the two different
568         ! body-fixed coordinate systems:
569 <      
569 >
570         xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
571         yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
572         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
515 –      
516 – #endif
573  
574 + #endif
575 +       xihat = xi / rij
576 +       yihat = yi / rij
577 +       zihat = zi / rij
578         xi2 = xi*xi
579         yi2 = yi*yi
580 <       zi2 = zi*zi
521 <      
522 <       proji = sqrt(xi2 + yi2)
523 <       proji3 = proji*proji*proji
524 <      
580 >       zi2 = zi*zi            
581         cti = zi / rij
582  
583 +       if (cti .gt. 1.0_dp) cti = 1.0_dp
584 +       if (cti .lt. -1.0_dp) cti = -1.0_dp
585 +
586         dctidx = - zi * xi / r3
587         dctidy = - zi * yi / r3
588         dctidz = 1.0d0 / rij - zi2 / r3
589 <       dctidux =  yi / rij
590 <       dctiduy = -xi / rij
591 <       dctiduz = 0.0d0
592 <      
589 >       dctidux = yi / rij ! - (zi * xi2) / r3
590 >       dctiduy = -xi / rij !- (zi * yi2) / r3
591 >       dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3
592 >
593 >       ! this is an attempt to try to truncate the singularity when
594 >       ! sin(theta) is near 0.0:
595 >
596 >       sti2 = 1.0_dp - cti*cti
597 >       if (dabs(sti2) .lt. 1.0d-12) then
598 >          proji = sqrt(rij * 1.0d-12)
599 >          dcpidx = 1.0d0 / proji
600 >          dcpidy = 0.0d0
601 >          dcpidux = xi / proji
602 >          dcpiduy = 0.0d0
603 >          dspidx = 0.0d0
604 >          dspidy = 1.0d0 / proji
605 >          dspidux = 0.0d0
606 >          dspiduy = yi / proji
607 >       else
608 >          proji = sqrt(xi2 + yi2)
609 >          proji3 = proji*proji*proji
610 >          dcpidx = 1.0d0 / proji - xi2 / proji3
611 >          dcpidy = - xi * yi / proji3
612 >          dcpidux = xi / proji - (xi2 * xi) / proji3
613 >          dcpiduy = - (xi * yi2) / proji3
614 >          dspidx = - xi * yi / proji3
615 >          dspidy = 1.0d0 / proji - yi2 / proji3
616 >          dspidux = - (yi * xi2) / proji3
617 >          dspiduy = yi / proji - (yi2 * yi) / proji3
618 >       endif
619 >
620         cpi = xi / proji
535 –       dcpidx = 1.0d0 / proji - xi2 / proji3
536 –       dcpidy = - xi * yi / proji3
621         dcpidz = 0.0d0
622 <       dcpidux = xi * yi * zi / proji3
623 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
540 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
541 <      
622 >       dcpiduz = 0.0d0
623 >
624         spi = yi / proji
543 –       dspidx = - xi * yi / proji3
544 –       dspidy = 1.0d0 / proji - yi2 / proji3
625         dspidz = 0.0d0
626 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
547 <       dspiduy = xi * yi * zi / proji3
548 <       dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
626 >       dspiduz = 0.0d0
627  
628         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
629              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 555 | Line 633 | contains  
633              CHEBYSHEV_TN, tm_i, dtm_i)
634         call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
635              CHEBYSHEV_UN, um_i, dum_i)
636 <      
636 >
637         sigma_i = 0.0d0
638         s_i = 0.0d0
639         eps_i = 0.0d0
# Line 589 | Line 667 | contains  
667               dPhuncdX = coeff * dtm_i(m) * dcpidx
668               dPhuncdY = coeff * dtm_i(m) * dcpidy
669               dPhuncdZ = coeff * dtm_i(m) * dcpidz
670 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
670 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
671               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
672               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
673            else
# Line 603 | Line 681 | contains  
681            endif
682  
683            sigma_i = sigma_i + plm_i(m,l)*Phunc
684 <
684 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux
685 > !!$          write(*,*) 'Phunc = ', Phunc
686            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
687                 Phunc * dlm_i(m,l) * dctidx
688            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
689                 Phunc * dlm_i(m,l) * dctidy
690            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
691                 Phunc * dlm_i(m,l) * dctidz
613 –          
692            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
693                 Phunc * dlm_i(m,l) * dctidux
694            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
695                 Phunc * dlm_i(m,l) * dctiduy
696            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
697                 Phunc * dlm_i(m,l) * dctiduz
698 <
698 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
699 > !!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
700 > !!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
701         end do
702  
703         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 625 | Line 705 | contains  
705            m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
706            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
707            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
708 <          
708 >
709            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
710               Phunc = coeff * tm_i(m)
711               dPhuncdX = coeff * dtm_i(m) * dcpidx
712               dPhuncdY = coeff * dtm_i(m) * dcpidy
713               dPhuncdZ = coeff * dtm_i(m) * dcpidz
714 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
714 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
715               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
716               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
717            else
# Line 645 | Line 725 | contains  
725            endif
726  
727            s_i = s_i + plm_i(m,l)*Phunc
728 <          
728 >
729            dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
730                 Phunc * dlm_i(m,l) * dctidx
731            dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
732                 Phunc * dlm_i(m,l) * dctidy
733            dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
734                 Phunc * dlm_i(m,l) * dctidz
735 <          
735 >
736            dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
737                 Phunc * dlm_i(m,l) * dctidux
738            dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 661 | Line 741 | contains  
741                 Phunc * dlm_i(m,l) * dctiduz      
742  
743         end do
744 <              
744 >
745         do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
746            l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
747            m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
748            coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
749            function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
750 <          
750 >
751            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
752               Phunc = coeff * tm_i(m)
753               dPhuncdX = coeff * dtm_i(m) * dcpidx
754               dPhuncdY = coeff * dtm_i(m) * dcpidy
755               dPhuncdZ = coeff * dtm_i(m) * dcpidz
756 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
756 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
757               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
758               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
759            else
# Line 687 | Line 767 | contains  
767            endif
768  
769            eps_i = eps_i + plm_i(m,l)*Phunc
770 <          
770 >
771            depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
772                 Phunc * dlm_i(m,l) * dctidx
773            depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
774                 Phunc * dlm_i(m,l) * dctidy
775            depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
776                 Phunc * dlm_i(m,l) * dctidz
777 <          
777 >
778            depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
779                 Phunc * dlm_i(m,l) * dctidux
780            depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 705 | Line 785 | contains  
785         end do
786  
787      endif
708 –      
709 –       ! now do j:
788  
789 +    ! now do j:
790 +
791      if (ShapeMap%Shapes(st2)%isLJ) then
792         sigma_j = ShapeMap%Shapes(st2)%sigma
793         s_j = ShapeMap%Shapes(st2)%sigma
# Line 731 | Line 811 | contains  
811         depsjduy = 0.0d0
812         depsjduz = 0.0d0
813      else
814 <      
814 >
815   #ifdef IS_MPI
816         ! rotate the inter-particle separation into the two different
817         ! body-fixed coordinate systems:
818         ! negative sign because this is the vector from j to i:
819 <      
819 >
820         xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
821         yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
822         zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
# Line 744 | Line 824 | contains  
824         ! rotate the inter-particle separation into the two different
825         ! body-fixed coordinate systems:
826         ! negative sign because this is the vector from j to i:
827 <      
827 >
828         xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
829         yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
830         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
831   #endif
832 <      
832 >
833 >       xjhat = xj / rij
834 >       yjhat = yj / rij
835 >       zjhat = zj / rij
836         xj2 = xj*xj
837         yj2 = yj*yj
838         zj2 = zj*zj
756 –      
757 –       projj = sqrt(xj2 + yj2)
758 –       projj3 = projj*projj*projj
759 –      
839         ctj = zj / rij
840 +
841 +       if (ctj .gt. 1.0_dp) ctj = 1.0_dp
842 +       if (ctj .lt. -1.0_dp) ctj = -1.0_dp
843 +
844         dctjdx = - zj * xj / r3
845         dctjdy = - zj * yj / r3
846         dctjdz = 1.0d0 / rij - zj2 / r3
847 <       dctjdux =  yj / rij
848 <       dctjduy = -xj / rij
849 <       dctjduz = 0.0d0
850 <      
847 >       dctjdux = yj / rij !- (zi * xj2) / r3
848 >       dctjduy = -xj / rij !- (zj * yj2) / r3
849 >       dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3
850 >
851 >       ! this is an attempt to try to truncate the singularity when
852 >       ! sin(theta) is near 0.0:
853 >
854 >       stj2 = 1.0_dp - ctj*ctj
855 >       if (dabs(stj2) .lt. 1.0d-12) then
856 >          projj = sqrt(rij * 1.0d-12)
857 >          dcpjdx = 1.0d0 / projj
858 >          dcpjdy = 0.0d0
859 >          dcpjdux = xj / projj
860 >          dcpjduy = 0.0d0
861 >          dspjdx = 0.0d0
862 >          dspjdy = 1.0d0 / projj
863 >          dspjdux = 0.0d0
864 >          dspjduy = yj / projj
865 >       else
866 >          projj = sqrt(xj2 + yj2)
867 >          projj3 = projj*projj*projj
868 >          dcpjdx = 1.0d0 / projj - xj2 / projj3
869 >          dcpjdy = - xj * yj / projj3
870 >          dcpjdux = xj / projj - (xj2 * xj) / projj3
871 >          dcpjduy = - (xj * yj2) / projj3
872 >          dspjdx = - xj * yj / projj3
873 >          dspjdy = 1.0d0 / projj - yj2 / projj3
874 >          dspjdux = - (yj * xj2) / projj3
875 >          dspjduy = yj / projj - (yj2 * yj) / projj3
876 >       endif
877 >
878         cpj = xj / projj
769 –       dcpjdx = 1.0d0 / projj - xj2 / projj3
770 –       dcpjdy = - xj * yj / projj3
879         dcpjdz = 0.0d0
880 <       dcpjdux = xj * yj * zj / projj3
881 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
774 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
775 <      
880 >       dcpjduz = 0.0d0
881 >
882         spj = yj / projj
777 –       dspjdx = - xj * yj / projj3
778 –       dspjdy = 1.0d0 / projj - yj2 / projj3
883         dspjdz = 0.0d0
884 <       dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
885 <       dspjduy = xj * yj * zj / projj3
886 <       dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
887 <      
884 >       dspjduz = 0.0d0
885 >
886 >
887 > !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
888 > !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
889         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
890              ShapeMap%Shapes(st2)%bigL, LMAX, &
891              plm_j, dlm_j)
892 <      
892 >
893         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
894              CHEBYSHEV_TN, tm_j, dtm_j)
895         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
896              CHEBYSHEV_UN, um_j, dum_j)
897 <      
897 >
898         sigma_j = 0.0d0
899         s_j = 0.0d0
900         eps_j = 0.0d0
# Line 823 | Line 928 | contains  
928               dPhuncdX = coeff * dtm_j(m) * dcpjdx
929               dPhuncdY = coeff * dtm_j(m) * dcpjdy
930               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
931 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
931 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
932               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
933               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
934            else
# Line 835 | Line 940 | contains  
940               dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
941               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
942            endif
943 <
943 >
944            sigma_j = sigma_j + plm_j(m,l)*Phunc
945 <          
945 >
946            dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
947                 Phunc * dlm_j(m,l) * dctjdx
948            dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
949                 Phunc * dlm_j(m,l) * dctjdy
950            dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
951                 Phunc * dlm_j(m,l) * dctjdz
952 <          
952 >
953            dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
954                 Phunc * dlm_j(m,l) * dctjdux
955            dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
# Line 865 | Line 970 | contains  
970               dPhuncdX = coeff * dtm_j(m) * dcpjdx
971               dPhuncdY = coeff * dtm_j(m) * dcpjdy
972               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
973 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
973 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
974               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
975               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
976            else
# Line 879 | Line 984 | contains  
984            endif
985  
986            s_j = s_j + plm_j(m,l)*Phunc
987 <          
987 >
988            dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
989                 Phunc * dlm_j(m,l) * dctjdx
990            dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
991                 Phunc * dlm_j(m,l) * dctjdy
992            dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
993                 Phunc * dlm_j(m,l) * dctjdz
994 <          
994 >
995            dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
996                 Phunc * dlm_j(m,l) * dctjdux
997            dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 920 | Line 1025 | contains  
1025               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1026            endif
1027  
1028 + !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1029 +
1030            eps_j = eps_j + plm_j(m,l)*Phunc
1031 <          
1031 >
1032            depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1033                 Phunc * dlm_j(m,l) * dctjdx
1034            depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1035                 Phunc * dlm_j(m,l) * dctjdy
1036            depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1037                 Phunc * dlm_j(m,l) * dctjdz
1038 <          
1038 >
1039            depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1040                 Phunc * dlm_j(m,l) * dctjdux
1041            depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 943 | Line 1050 | contains  
1050      ! phew, now let's assemble the potential energy:
1051  
1052      sigma = 0.5*(sigma_i + sigma_j)
1053 <
1053 > !    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1054      dsigmadxi = 0.5*dsigmaidx
1055      dsigmadyi = 0.5*dsigmaidy
1056      dsigmadzi = 0.5*dsigmaidz
# Line 975 | Line 1082 | contains  
1082      dsduzj = 0.5*dsjduz
1083  
1084      eps = sqrt(eps_i * eps_j)
1085 <
1085 > !!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1086 > !!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1087 > !!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1088      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1089      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1090      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 989 | Line 1098 | contains  
1098      depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1099      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1100      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1101 <    
1101 >
1102 > !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1103 >
1104 > !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1105 > !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1106 > !!$
1107 > !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1108 >
1109      rtdenom = rij-sigma+s
1110      rt = s / rtdenom
1111  
1112 <    drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1113 <    drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1114 <    drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1115 <    drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1116 <    drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1117 <    drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1118 <    drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1119 <    drtdyj = (dsdyj + rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1120 <    drtdzj = (dsdzj + rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1121 <    drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1122 <    drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1123 <    drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1124 <    
1112 >    drtdxi = (dsdxi - rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1113 >    drtdyi = (dsdyi - rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1114 >    drtdzi = (dsdzi - rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1115 >    drtduxi = (dsduxi - rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1116 >    drtduyi = (dsduyi - rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1117 >    drtduzi = (dsduzi - rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1118 >    drtdxj = (dsdxj - rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1119 >    drtdyj = (dsdyj - rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1120 >    drtdzj = (dsdzj - rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1121 >    drtduxj = (dsduxj - rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1122 >    drtduyj = (dsduyj - rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1123 >    drtduzj = (dsduzj - rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1124 >
1125 > !!$    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1126 > !!$    write(*,*) 'drtdu_j = ', drtduxj, drtduyj, drtduzj
1127 >
1128      rt2 = rt*rt
1129      rt3 = rt2*rt
1130      rt5 = rt2*rt3
# Line 1014 | Line 1133 | contains  
1133      rt12 = rt6*rt6
1134      rt126 = rt12 - rt6
1135  
1136 +    pot_temp = 4.0d0 * eps * rt126
1137 +
1138 +    vpair = vpair + pot_temp
1139      if (do_pot) then
1140   #ifdef IS_MPI
1141 <       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1142 <       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1141 >       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1142 >       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1143   #else
1144 <       pot = pot + 4.0d0*eps*rt126*sw
1144 >       pot = pot + pot_temp*sw
1145   #endif
1146      endif
1147 <    
1147 >
1148 > !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1149 >
1150      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1151      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1152      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1036 | Line 1160 | contains  
1160      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1161      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1162      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1163 <
1163 > !!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1164      ! do the torques first since they are easy:
1165      ! remember that these are still in the body fixed axes
1166  
1167 <    txi = dvduxi * sw
1168 <    tyi = dvduyi * sw
1169 <    tzi = dvduzi * sw
1046 <
1047 <    txj = dvduxj * sw
1048 <    tyj = dvduyj * sw
1049 <    tzj = dvduzj * sw
1167 >    txi = 0.0d0
1168 >    tyi = 0.0d0
1169 >    tzi = 0.0d0
1170  
1171 +    txj = 0.0d0
1172 +    tyj = 0.0d0
1173 +    tzj = 0.0d0
1174 +
1175 +    txi = (dvduyi - dvduzi) * sw
1176 +    tyi = (dvduzi - dvduxi) * sw
1177 +    tzi = (dvduxi - dvduyi) * sw
1178 +
1179 +    txj = (dvduyj - dvduzj) * sw
1180 +    tyj = (dvduzj - dvduxj) * sw
1181 +    tzj = (dvduxj - dvduyj) * sw
1182 +
1183 + !!$    txi = dvduxi * sw
1184 + !!$    tyi = dvduyi * sw
1185 + !!$    tzi = dvduzi * sw
1186 + !!$
1187 + !!$    txj = dvduxj * sw
1188 + !!$    tyj = dvduyj * sw
1189 + !!$    tzj = dvduzj * sw
1190 +
1191 +    write(*,*) 't1 = ', txi, tyi, tzi
1192 +    write(*,*) 't2 = ', txj, tyj, tzj
1193 +
1194      ! go back to lab frame using transpose of rotation matrix:
1195 <    
1195 >
1196   #ifdef IS_MPI
1197      t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
1198           a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
# Line 1057 | Line 1200 | contains  
1200           a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
1201      t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
1202           a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
1203 <    
1203 >
1204      t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
1205           a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
1206      t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
1207 <            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1207 >         a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1208      t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
1209           a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
1210   #else
1211      t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
1212      t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
1213      t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
1214 <    
1214 >
1215      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1216      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1217      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1218 +
1219   #endif
1220      ! Now, on to the forces:
1221 <    
1221 >
1222      ! first rotate the i terms back into the lab frame:
1079 –    
1080 –    fxi = dvdxi * sw
1081 –    fyi = dvdyi * sw
1082 –    fzi = dvdzi * sw
1223  
1224 <    fxj = dvdxj * sw
1225 <    fyj = dvdyj * sw
1226 <    fzj = dvdzj * sw
1224 >    fxi = -dvdxi * sw
1225 >    fyi = -dvdyi * sw
1226 >    fzi = -dvdzi * sw
1227  
1228 +    fxj = -dvdxj * sw
1229 +    fyj = -dvdyj * sw
1230 +    fzj = -dvdzj * sw
1231 +
1232 +
1233   #ifdef IS_MPI
1234      fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1235      fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
# Line 1097 | Line 1242 | contains  
1242      fxii = a(1,atom1)*fxi + a(4,atom1)*fyi + a(7,atom1)*fzi
1243      fyii = a(2,atom1)*fxi + a(5,atom1)*fyi + a(8,atom1)*fzi
1244      fzii = a(3,atom1)*fxi + a(6,atom1)*fyi + a(9,atom1)*fzi
1245 <    
1245 >
1246      fxjj = a(1,atom2)*fxj + a(4,atom2)*fyj + a(7,atom2)*fzj
1247      fyjj = a(2,atom2)*fxj + a(5,atom2)*fyj + a(8,atom2)*fzj
1248      fzjj = a(3,atom2)*fxj + a(6,atom2)*fyj + a(9,atom2)*fzj
# Line 1106 | Line 1251 | contains  
1251      fxij = -fxii
1252      fyij = -fyii
1253      fzij = -fzii
1254 <    
1254 >
1255      fxji = -fxjj
1256      fyji = -fyjj
1257      fzji = -fzjj
1258  
1259 <    fxradial = fxii + fxji
1260 <    fyradial = fyii + fyji
1261 <    fzradial = fzii + fzji
1262 <
1259 >    fxradial = (fxii + fxji)
1260 >    fyradial = (fyii + fyji)
1261 >    fzradial = (fzii + fzji)
1262 > !!$    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1263   #ifdef IS_MPI
1264      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1265      f_Row(2,atom1) = f_Row(2,atom1) + fyradial
1266      f_Row(3,atom1) = f_Row(3,atom1) + fzradial
1267 <    
1267 >
1268      f_Col(1,atom2) = f_Col(1,atom2) - fxradial
1269      f_Col(2,atom2) = f_Col(2,atom2) - fyradial
1270      f_Col(3,atom2) = f_Col(3,atom2) - fzradial
# Line 1127 | Line 1272 | contains  
1272      f(1,atom1) = f(1,atom1) + fxradial
1273      f(2,atom1) = f(2,atom1) + fyradial
1274      f(3,atom1) = f(3,atom1) + fzradial
1275 <    
1275 >
1276      f(1,atom2) = f(1,atom2) - fxradial
1277      f(2,atom2) = f(2,atom2) - fyradial
1278      f(3,atom2) = f(3,atom2) - fzradial
# Line 1140 | Line 1285 | contains  
1285      id1 = atom1
1286      id2 = atom2
1287   #endif
1288 <    
1288 >
1289      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1290 <      
1290 >
1291         fpair(1) = fpair(1) + fxradial
1292         fpair(2) = fpair(2) + fyradial
1293         fpair(3) = fpair(3) + fzradial
1294 <      
1294 >
1295      endif
1296 <    
1296 >
1297    end subroutine do_shape_pair
1298 <    
1298 >
1299    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1300  
1301      ! Purpose: Compute the associated Legendre functions
# Line 1168 | Line 1313 | contains  
1313      !
1314      ! The original Fortran77 codes can be found here:
1315      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1316 <    
1316 >
1317      real (kind=dp), intent(in) :: x
1318      integer, intent(in) :: l, m, lmax
1319      real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
# Line 1185 | Line 1330 | contains  
1330  
1331      ! start with 0,0:
1332      PLM(0,0) = 1.0D0
1333 <  
1333 >
1334      ! x = +/- 1 functions are easy:
1335      IF (abs(X).EQ.1.0D0) THEN
1336         DO I = 1, m
# Line 1240 | Line 1385 | contains  
1385  
1386  
1387    subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1388 <  
1388 >
1389      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1390      !          or Ln(x) or Hn(x), and their derivatives
1391      ! Input :  function_type --- Function code
# Line 1259 | Line 1404 | contains  
1404      !
1405      ! The original Fortran77 codes can be found here:
1406      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1407 <  
1407 >
1408      real(kind=8), intent(in) :: x
1409      integer, intent(in):: m, mmax
1410      integer, intent(in):: function_type
1411      real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1412 <  
1412 >
1413      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1414      integer :: k
1415  
# Line 1310 | Line 1455 | contains  
1455  
1456  
1457      RETURN
1458 <    
1458 >
1459    end subroutine Orthogonal_Polynomial
1315 –  
1316 – end module shapes
1460  
1461 < subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1462 <     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1320 <     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1321 <     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1322 <     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1323 <     myATID, status)
1461 >  subroutine deallocateShapes(this)
1462 >    type(Shape), pointer :: this
1463  
1464 <  use definitions
1465 <  use shapes, only: newShapeType
1466 <  
1467 <  integer :: nContactFuncs
1329 <  integer :: nRangeFuncs
1330 <  integer :: nStrengthFuncs
1331 <  integer :: status
1332 <  integer :: myATID
1333 <  
1334 <  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1335 <  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1336 <  integer, dimension(nContactFuncs) :: ContactFunctionType        
1337 <  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1338 <  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1339 <  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1340 <  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1341 <  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1342 <  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1343 <  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1344 <  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1345 <  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1346 <  
1347 <  call newShapeType(nContactFuncs, ContactFuncLValue, &
1348 <       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1349 <       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1350 <       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1351 <       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1352 <       myATID, status)
1464 >    if (associated( this%ContactFuncLValue)) then
1465 >       deallocate(this%ContactFuncLValue)
1466 >       this%ContactFuncLValue => null()
1467 >    end if
1468  
1469 <  return
1470 < end subroutine makeShape
1469 >    if (associated( this%ContactFuncMValue)) then
1470 >       deallocate( this%ContactFuncMValue)
1471 >       this%ContactFuncMValue => null()
1472 >    end if
1473 >    if (associated( this%ContactFunctionType)) then
1474 >       deallocate(this%ContactFunctionType)
1475 >       this%ContactFunctionType => null()
1476 >    end if
1477  
1478 < subroutine completeShapeFF(status)
1478 >    if (associated( this%ContactFuncCoefficient)) then
1479 >       deallocate(this%ContactFuncCoefficient)
1480 >       this%ContactFuncCoefficient => null()
1481 >    end if
1482  
1483 <  use shapes, only: complete_Shape_FF
1483 >    if (associated( this%RangeFuncLValue)) then
1484 >       deallocate(this%RangeFuncLValue)
1485 >       this%RangeFuncLValue => null()
1486 >    end if
1487 >    if (associated( this%RangeFuncMValue)) then
1488 >       deallocate( this%RangeFuncMValue)
1489 >       this%RangeFuncMValue => null()
1490 >    end if
1491  
1492 <  integer, intent(out)  :: status
1493 <  integer :: myStatus
1492 >    if (associated( this%RangeFunctionType)) then
1493 >       deallocate( this%RangeFunctionType)
1494 >       this%RangeFunctionType => null()
1495 >    end if
1496 >    if (associated( this%RangeFuncCoefficient)) then
1497 >       deallocate(this%RangeFuncCoefficient)
1498 >       this%RangeFuncCoefficient => null()
1499 >    end if
1500  
1501 <  myStatus = 0
1501 >    if (associated( this%StrengthFuncLValue)) then
1502 >       deallocate(this%StrengthFuncLValue)
1503 >       this%StrengthFuncLValue => null()
1504 >    end if
1505  
1506 <  call complete_Shape_FF(myStatus)
1506 >    if (associated( this%StrengthFuncMValue )) then
1507 >       deallocate(this%StrengthFuncMValue)
1508 >       this%StrengthFuncMValue => null()
1509 >    end if
1510  
1511 <  status = myStatus
1511 >    if(associated( this%StrengthFunctionType)) then
1512 >       deallocate(this%StrengthFunctionType)
1513 >       this%StrengthFunctionType => null()
1514 >    end if
1515 >    if (associated( this%StrengthFuncCoefficient )) then
1516 >       deallocate(this%StrengthFuncCoefficient)
1517 >       this%StrengthFuncCoefficient => null()
1518 >    end if
1519 >  end subroutine deallocateShapes
1520  
1521 <  return
1522 < end subroutine completeShapeFF
1521 >  subroutine destroyShapeTypes
1522 >    integer :: i
1523 >    type(Shape), pointer :: thisShape
1524  
1525 +    ! First walk through and kill the shape
1526 +    do i = 1,ShapeMap%n_shapes
1527 +       thisShape => ShapeMap%Shapes(i)
1528 +       call deallocateShapes(thisShape)
1529 +    end do
1530 +
1531 +    ! set shape map to starting values
1532 +    ShapeMap%n_shapes = 0
1533 +    ShapeMap%currentShape = 0
1534 +
1535 +    if (associated(ShapeMap%Shapes)) then
1536 +       deallocate(ShapeMap%Shapes)
1537 +       ShapeMap%Shapes => null()
1538 +    end if
1539 +
1540 +    if (associated(ShapeMap%atidToShape)) then
1541 +       deallocate(ShapeMap%atidToShape)
1542 +       ShapeMap%atidToShape => null()
1543 +    end if
1544 +
1545 +
1546 +  end subroutine destroyShapeTypes
1547 +
1548 +
1549 + end module shapes

Diff Legend

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