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 1633 by gezelter, Fri Oct 22 20:22:48 2004 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 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 25 | Line 67 | module shapes
67  
68    public :: do_shape_pair
69    public :: newShapeType
70 +  public :: complete_Shape_FF
71 +  public :: destroyShapeTypes
72  
29 –
73    type, private :: Shape
74       integer :: atid
75       integer :: nContactFuncs
# Line 50 | 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()
102    end type ShapeList
103 <  
103 >
104    type(ShapeList), save :: ShapeMap
105  
106    integer :: lmax
64 –  real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
65 –  real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
66 –  real (kind=dp), allocatable, dimension(:) :: tm_j, dtm_j, um_j, dum_j
107  
108   contains  
109  
# Line 72 | Line 112 | contains  
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115 <       myAtid, status)
115 >       myATID, status)
116  
117      integer :: nContactFuncs
118      integer :: nRangeFuncs
119      integer :: nStrengthFuncs
120      integer :: shape_ident
121      integer :: status
122 <    integer :: myAtid
122 >    integer :: myATID
123      integer :: bigL
124      integer :: bigM
125      integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
# Line 104 | Line 144 | contains  
144  
145         call getMatchingElementList(atypes, "is_Shape", .true., &
146              nShapeTypes, MatchList)
147 <      
147 >
148         call getMatchingElementList(atypes, "is_LennardJones", .true., &
149              nLJTypes, MatchList)
150 <      
150 >
151         ShapeMap%n_shapes = nShapeTypes + nLJTypes
152 <      
152 >
153         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
154 <      
154 >
155         ntypes = getSize(atypes)
156 <      
157 <       allocate(ShapeMap%atidToShape(ntypes))
156 >
157 >       allocate(ShapeMap%atidToShape(0:ntypes))
158      end if
159 <    
159 >
160      ShapeMap%currentShape = ShapeMap%currentShape + 1
161      current = ShapeMap%currentShape
162  
# Line 127 | Line 167 | contains  
167         return
168      endif
169  
170 <    call getElementProperty(atypes, myAtid, "c_ident", me)
170 >    call getElementProperty(atypes, myATID, 'c_ident', me)
171 >
172      ShapeMap%atidToShape(me)                         = current
173      ShapeMap%Shapes(current)%atid                    = me
174      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
# Line 148 | Line 189 | contains  
189  
190      bigL = -1
191      bigM = -1
192 <    
192 >
193      do j = 1, ShapeMap%Shapes(current)%nContactFuncs
194         if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
195            bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
# Line 186 | Line 227 | contains  
227      type(Shape), intent(inout) :: myShape
228      integer, intent(out) :: stat
229      integer :: alloc_stat
230 <
230 >
231 >    stat = 0
232      if (associated(myShape%contactFuncLValue)) then
233         deallocate(myShape%contactFuncLValue)
234      endif
# Line 252 | Line 294 | contains  
294         stat = -1
295         return
296      endif
297 <    
297 >
298      if (associated(myShape%strengthFuncLValue)) then
299         deallocate(myShape%strengthFuncLValue)
300      endif
# Line 286 | Line 328 | contains  
328         return
329      endif
330  
331 +    return
332 +
333    end subroutine allocateShape
334 <    
335 <  subroutine init_Shape_FF(status)
334 >
335 >  subroutine complete_Shape_FF(status)
336      integer :: status
337      integer :: i, j, l, m, lm, function_type
338 <    real(kind=dp) :: bigSigma, myBigSigma, thisSigma, coeff, Phunc, spi
295 <    real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP, sigma
338 >    real(kind=dp) :: thisDP, sigma
339      integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
340      logical :: thisProperty
341  
299 –    Pi = 4.0d0 * datan(1.0d0)
300 –
342      status = 0
343      if (ShapeMap%currentShape == 0) then
344         call handleError("init_Shape_FF", "No members in ShapeMap")
345         status = -1
346         return
347      end if
307 –
308 –    bigSigma = 0.0d0
309 –    do i = 1, ShapeMap%currentShape
310 –
311 –       ! Scan over theta and phi to
312 –       ! find the largest contact in any direction....
313 –
314 –       myBigSigma = 0.0d0
315 –
316 –       do iTheta = 0, nSteps
317 –          theta = (Pi/2.0d0)*(dble(iTheta)/dble(nSteps))
318 –          costheta = cos(theta)
319 –
320 –          call Associated_Legendre(costheta, ShapeMap%Shapes(i)%bigL, &
321 –               ShapeMap%Shapes(i)%bigM, lmax, plm_i, dlm_i)
322 –                  
323 –          do iPhi = 0, nSteps
324 –             phi = -Pi  + 2.0d0 * Pi * (dble(iPhi)/dble(nSteps))
325 –             cpi = cos(phi)
326 –             spi = sin(phi)
327 –            
328 –             call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
329 –                  CHEBYSHEV_TN, tm_i, dtm_i)
330 –             call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
331 –                  CHEBYSHEV_UN, um_i, dum_i)
348  
333 –             thisSigma = 0.0d0
334 –
335 –             do lm = 1, ShapeMap%Shapes(i)%nContactFuncs                
336 –
337 –                l = ShapeMap%Shapes(i)%ContactFuncLValue(lm)
338 –                m = ShapeMap%Shapes(i)%ContactFuncMValue(lm)
339 –                coeff = ShapeMap%Shapes(i)%ContactFuncCoefficient(lm)
340 –                function_type = ShapeMap%Shapes(i)%ContactFunctionType(lm)
341 –
342 –                if ((function_type .eq. SH_COS).or.(m.eq.0)) then
343 –                   Phunc = coeff * tm_i(m)
344 –                else
345 –                   Phunc = coeff * spi * um_i(m-1)
346 –                endif
347 –                
348 –                thisSigma = thisSigma + plm_i(l,m)*Phunc
349 –             enddo
350 –
351 –             if (thisSigma.gt.myBigSigma) myBigSigma = thisSigma
352 –          enddo
353 –       enddo
354 –
355 –       if (myBigSigma.gt.bigSigma) bigSigma = myBigSigma
356 –    enddo
357 –
349      nAtypes = getSize(atypes)
350  
351      if (nAtypes == 0) then
# Line 362 | Line 353 | contains  
353         return
354      end if
355  
356 <    do i = 1, nAtypes
357 <      
356 >    ! atypes comes from c side
357 >    do i = 0, nAtypes
358 >
359         call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
360 <      
360 >
361         if (thisProperty) then
362 <          
362 >
363            ShapeMap%currentShape = ShapeMap%currentShape + 1
364            current = ShapeMap%currentShape
365  
# Line 378 | Line 370 | contains  
370            ShapeMap%Shapes(current)%isLJ = .true.
371  
372            ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
373 <          sigma = getSigma(thisIP)
374 <          ShapeMap%Shapes(current)%sigma = sigma
383 <          if (sigma .gt. bigSigma) bigSigma = thisDP
384 <          
373 >          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
374 >
375         endif
376 <      
376 >
377      end do
378  
379      haveShapeMap = .true.
380 <    
381 <  end subroutine init_Shape_FF
382 <    
380 >
381 >  end subroutine complete_Shape_FF
382 >
383    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
384         pot, A, f, t, do_pot)
385 <    
385 >
386 >    INTEGER, PARAMETER:: LMAX         = 64
387 >    INTEGER, PARAMETER:: MMAX         = 64
388 >
389      integer, intent(in) :: atom1, atom2
390      real (kind=dp), intent(inout) :: rij, r2
391      real (kind=dp), dimension(3), intent(in) :: d
392      real (kind=dp), dimension(3), intent(inout) :: fpair
393 <    real (kind=dp) :: pot, vpair, sw
393 >    real (kind=dp) :: pot, vpair, sw, dswdr
394      real (kind=dp), dimension(9,nLocal) :: A
395      real (kind=dp), dimension(3,nLocal) :: f
396      real (kind=dp), dimension(3,nLocal) :: t
# Line 408 | Line 401 | contains  
401      integer :: l, m, lm, id1, id2, localError, function_type
402      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
403      real (kind=dp) :: coeff
404 +    real (kind=dp) :: pot_temp
405  
406      real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
407      real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
# Line 426 | Line 420 | contains  
420  
421      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
422  
423 +    real (kind=dp) :: sti2, stj2
424 +
425      real (kind=dp) :: proji, proji3, projj, projj3
426      real (kind=dp) :: cti, ctj, cpi, cpj, spi, spj
427      real (kind=dp) :: Phunc, sigma, s, eps, rtdenom, rt
# Line 457 | Line 453 | contains  
453      real (kind=dp) :: dsduxi, dsduyi, dsduzi
454      real (kind=dp) :: dsdxj, dsdyj, dsdzj
455      real (kind=dp) :: dsduxj, dsduyj, dsduzj
456 <    
456 >
457      real (kind=dp) :: depsdxi, depsdyi, depsdzi
458      real (kind=dp) :: depsduxi, depsduyi, depsduzi
459      real (kind=dp) :: depsdxj, depsdyj, depsdzj
# Line 484 | Line 480 | contains  
480      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
481      real (kind=dp) :: fxradial, fyradial, fzradial
482  
483 +    real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
484 +    real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
485 +    real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
486 +    real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
487 +
488      if (.not.haveShapeMap) then
489         call handleError("calc_shape", "NO SHAPEMAP!!!!")
490         return      
491      endif
492 <    
492 >
493      !! We assume that the rotation matrices have already been calculated
494      !! and placed in the A array.
495 <        
495 >
496      r3 = r2*rij
497      r5 = r3*r2
498 <    
498 >
499      drdxi = -d(1) / rij
500      drdyi = -d(2) / rij
501      drdzi = -d(3) / rij
# Line 502 | Line 503 | contains  
503      drdxj = d(1) / rij
504      drdyj = d(2) / rij
505      drdzj = d(3) / rij
506 <    
506 >
507      ! find the atom type id (atid) for each atom:
508   #ifdef IS_MPI
509      atid1 = atid_Row(atom1)
# Line 513 | Line 514 | contains  
514   #endif
515  
516      ! use the atid to find the shape type (st) for each atom:
516 –
517      st1 = ShapeMap%atidToShape(atid1)
518      st2 = ShapeMap%atidToShape(atid2)
519 <    
519 >
520      if (ShapeMap%Shapes(st1)%isLJ) then
521 +
522         sigma_i = ShapeMap%Shapes(st1)%sigma
523         s_i = ShapeMap%Shapes(st1)%sigma
524         eps_i = ShapeMap%Shapes(st1)%epsilon
# Line 544 | Line 545 | contains  
545   #ifdef IS_MPI
546         ! rotate the inter-particle separation into the two different
547         ! body-fixed coordinate systems:
548 <      
548 >
549         xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
550         yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
551         zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
552 <      
552 >
553   #else
554         ! rotate the inter-particle separation into the two different
555         ! body-fixed coordinate systems:
556 <      
556 >
557         xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
558         yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
559         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
560 <      
560 >
561   #endif
562 <      
562 >
563         xi2 = xi*xi
564         yi2 = yi*yi
565 <       zi2 = zi*zi
565 <      
566 <       proji = sqrt(xi2 + yi2)
567 <       proji3 = proji*proji*proji
568 <      
565 >       zi2 = zi*zi            
566         cti = zi / rij
567 +
568 +       if (cti .gt. 1.0_dp) cti = 1.0_dp
569 +       if (cti .lt. -1.0_dp) cti = -1.0_dp
570 +
571         dctidx = - zi * xi / r3
572         dctidy = - zi * yi / r3
573         dctidz = 1.0d0 / rij - zi2 / r3
574 <       dctidux =  yi / rij
575 <       dctiduy = -xi / rij
576 <       dctiduz = 0.0d0
577 <      
574 >       dctidux = - (zi * xi2) / r3
575 >       dctiduy = - (zi * yi2) / r3
576 >       dctiduz = zi / rij - (zi2 * zi) / r3
577 >
578 >       ! this is an attempt to try to truncate the singularity when
579 >       ! sin(theta) is near 0.0:
580 >
581 >       sti2 = 1.0_dp - cti*cti
582 >       if (dabs(sti2) .lt. 1.0d-12) then
583 >          proji = sqrt(rij * 1.0d-12)
584 >          dcpidx = 1.0d0 / proji
585 >          dcpidy = 0.0d0
586 >          dcpidux = xi / proji
587 >          dcpiduy = 0.0d0
588 >          dspidx = 0.0d0
589 >          dspidy = 1.0d0 / proji
590 >          dspidux = 0.0d0
591 >          dspiduy = yi / proji
592 >       else
593 >          proji = sqrt(xi2 + yi2)
594 >          proji3 = proji*proji*proji
595 >          dcpidx = 1.0d0 / proji - xi2 / proji3
596 >          dcpidy = - xi * yi / proji3
597 >          dcpidux = xi / proji - (xi2 * xi) / proji3
598 >          dcpiduy = - (xi * yi2) / proji3
599 >          dspidx = - xi * yi / proji3
600 >          dspidy = 1.0d0 / proji - yi2 / proji3
601 >          dspidux = - (yi * xi2) / proji3
602 >          dspiduy = yi / proji - (yi2 * yi) / proji3
603 >       endif
604 >
605         cpi = xi / proji
578 –       dcpidx = 1.0d0 / proji - xi2 / proji3
579 –       dcpidy = - xi * yi / proji3
606         dcpidz = 0.0d0
607 <       dcpidux = xi * yi * zi / proji3
608 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
583 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
584 <      
607 >       dcpiduz = 0.0d0
608 >
609         spi = yi / proji
586 –       dspidx = - xi * yi / proji3
587 –       dspidy = 1.0d0 / proji - yi2 / proji3
610         dspidz = 0.0d0
611 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
590 <       dspiduy = xi * yi * zi / proji3
591 <       dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
611 >       dspiduz = 0.0d0
612  
613 <       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
614 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
613 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
614 >            ShapeMap%Shapes(st1)%bigL, LMAX, &
615 >            plm_i, dlm_i)
616  
617 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
617 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
618              CHEBYSHEV_TN, tm_i, dtm_i)
619 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
619 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
620              CHEBYSHEV_UN, um_i, dum_i)
621 <      
621 >
622         sigma_i = 0.0d0
623         s_i = 0.0d0
624         eps_i = 0.0d0
# Line 644 | Line 665 | contains  
665               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
666            endif
667  
668 <          sigma_i = sigma_i + plm_i(l,m)*Phunc
648 <          
649 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
650 <               Phunc * dlm_i(l,m) * dctidx
651 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
652 <               Phunc * dlm_i(l,m) * dctidy
653 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
654 <               Phunc * dlm_i(l,m) * dctidz
655 <          
656 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
657 <               Phunc * dlm_i(l,m) * dctidux
658 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
659 <               Phunc * dlm_i(l,m) * dctiduy
660 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
661 <               Phunc * dlm_i(l,m) * dctiduz
668 >          sigma_i = sigma_i + plm_i(m,l)*Phunc
669  
670 +          dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
671 +               Phunc * dlm_i(m,l) * dctidx
672 +          dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
673 +               Phunc * dlm_i(m,l) * dctidy
674 +          dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
675 +               Phunc * dlm_i(m,l) * dctidz
676 +
677 +          dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
678 +               Phunc * dlm_i(m,l) * dctidux
679 +          dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
680 +               Phunc * dlm_i(m,l) * dctiduy
681 +          dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
682 +               Phunc * dlm_i(m,l) * dctiduz
683 +
684         end do
685  
686         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 667 | Line 688 | contains  
688            m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
689            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
690            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
691 <          
691 >
692            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
693               Phunc = coeff * tm_i(m)
694               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 686 | Line 707 | contains  
707               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
708            endif
709  
710 <          s_i = s_i + plm_i(l,m)*Phunc
690 <          
691 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
692 <               Phunc * dlm_i(l,m) * dctidx
693 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
694 <               Phunc * dlm_i(l,m) * dctidy
695 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
696 <               Phunc * dlm_i(l,m) * dctidz
697 <          
698 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
699 <               Phunc * dlm_i(l,m) * dctidux
700 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
701 <               Phunc * dlm_i(l,m) * dctiduy
702 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
703 <               Phunc * dlm_i(l,m) * dctiduz      
710 >          s_i = s_i + plm_i(m,l)*Phunc
711  
712 +          dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
713 +               Phunc * dlm_i(m,l) * dctidx
714 +          dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
715 +               Phunc * dlm_i(m,l) * dctidy
716 +          dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
717 +               Phunc * dlm_i(m,l) * dctidz
718 +
719 +          dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
720 +               Phunc * dlm_i(m,l) * dctidux
721 +          dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
722 +               Phunc * dlm_i(m,l) * dctiduy
723 +          dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
724 +               Phunc * dlm_i(m,l) * dctiduz      
725 +
726         end do
727 <              
727 >
728         do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
729            l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
730            m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
731            coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
732            function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
733 <          
733 >
734            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
735               Phunc = coeff * tm_i(m)
736               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 728 | Line 749 | contains  
749               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
750            endif
751  
752 <          eps_i = eps_i + plm_i(l,m)*Phunc
732 <          
733 <          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
734 <               Phunc * dlm_i(l,m) * dctidx
735 <          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
736 <               Phunc * dlm_i(l,m) * dctidy
737 <          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
738 <               Phunc * dlm_i(l,m) * dctidz
739 <          
740 <          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
741 <               Phunc * dlm_i(l,m) * dctidux
742 <          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
743 <               Phunc * dlm_i(l,m) * dctiduy
744 <          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
745 <               Phunc * dlm_i(l,m) * dctiduz      
752 >          eps_i = eps_i + plm_i(m,l)*Phunc
753  
754 +          depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
755 +               Phunc * dlm_i(m,l) * dctidx
756 +          depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
757 +               Phunc * dlm_i(m,l) * dctidy
758 +          depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
759 +               Phunc * dlm_i(m,l) * dctidz
760 +
761 +          depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
762 +               Phunc * dlm_i(m,l) * dctidux
763 +          depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
764 +               Phunc * dlm_i(m,l) * dctiduy
765 +          depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
766 +               Phunc * dlm_i(m,l) * dctiduz      
767 +
768         end do
769  
770      endif
750 –      
751 –       ! now do j:
771  
772 +    ! now do j:
773 +
774      if (ShapeMap%Shapes(st2)%isLJ) then
775         sigma_j = ShapeMap%Shapes(st2)%sigma
776         s_j = ShapeMap%Shapes(st2)%sigma
# Line 773 | Line 794 | contains  
794         depsjduy = 0.0d0
795         depsjduz = 0.0d0
796      else
797 <      
797 >
798   #ifdef IS_MPI
799         ! rotate the inter-particle separation into the two different
800         ! body-fixed coordinate systems:
801         ! negative sign because this is the vector from j to i:
802 <      
802 >
803         xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
804         yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
805         zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
# Line 786 | Line 807 | contains  
807         ! rotate the inter-particle separation into the two different
808         ! body-fixed coordinate systems:
809         ! negative sign because this is the vector from j to i:
810 <      
810 >
811         xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
812         yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
813         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
814   #endif
815 <      
815 >
816         xj2 = xj*xj
817         yj2 = yj*yj
818         zj2 = zj*zj
798 –      
799 –       projj = sqrt(xj2 + yj2)
800 –       projj3 = projj*projj*projj
801 –      
819         ctj = zj / rij
820 +
821 +       if (ctj .gt. 1.0_dp) ctj = 1.0_dp
822 +       if (ctj .lt. -1.0_dp) ctj = -1.0_dp
823 +
824         dctjdx = - zj * xj / r3
825         dctjdy = - zj * yj / r3
826         dctjdz = 1.0d0 / rij - zj2 / r3
827 <       dctjdux =  yj / rij
828 <       dctjduy = -xj / rij
829 <       dctjduz = 0.0d0
830 <      
831 <       cpj = xj / projj
832 <       dcpjdx = 1.0d0 / projj - xj2 / projj3
833 <       dcpjdy = - xj * yj / projj3
834 <       dcpjdz = 0.0d0
835 <       dcpjdux = xj * yj * zj / projj3
836 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
837 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
838 <      
839 <       spj = yj / projj
840 <       dspjdx = - xj * yj / projj3
841 <       dspjdy = 1.0d0 / projj - yj2 / projj3
842 <       dspjdz = 0.0d0
843 <       dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
844 <       dspjduy = xj * yj * zj / projj3
845 <       dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
846 <      
847 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
848 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
849 <      
850 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
851 <            CHEBYSHEV_TN, tm_j, dtm_j)
852 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
827 >       dctjdux = - (zi * xj2) / r3
828 >       dctjduy = - (zj * yj2) / r3
829 >       dctjduz = zj / rij - (zj2 * zj) / r3
830 >
831 >       ! this is an attempt to try to truncate the singularity when
832 >       ! sin(theta) is near 0.0:
833 >
834 >       stj2 = 1.0_dp - ctj*ctj
835 >       if (dabs(stj2) .lt. 1.0d-12) then
836 >          projj = sqrt(rij * 1.0d-12)
837 >          dcpjdx = 1.0d0 / projj
838 >          dcpjdy = 0.0d0
839 >          dcpjdux = xj / projj
840 >          dcpjduy = 0.0d0
841 >          dspjdx = 0.0d0
842 >          dspjdy = 1.0d0 / projj
843 >          dspjdux = 0.0d0
844 >          dspjduy = yj / projj
845 >       else
846 >          projj = sqrt(xj2 + yj2)
847 >          projj3 = projj*projj*projj
848 >          dcpjdx = 1.0d0 / projj - xj2 / projj3
849 >          dcpjdy = - xj * yj / projj3
850 >          dcpjdux = xj / projj - (xj2 * xj) / projj3
851 >          dcpjduy = - (xj * yj2) / projj3
852 >          dspjdx = - xj * yj / projj3
853 >          dspjdy = 1.0d0 / projj - yj2 / projj3
854 >          dspjdux = - (yj * xj2) / projj3
855 >          dspjduy = yj / projj - (yj2 * yj) / projj3
856 >       endif
857 >
858 >       cpj = xj / projj
859 >       dcpjdz = 0.0d0
860 >       dcpjduz = 0.0d0
861 >
862 >       spj = yj / projj
863 >       dspjdz = 0.0d0
864 >       dspjduz = 0.0d0
865 >
866 >
867 >       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
868 >       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
869 >       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
870 >            ShapeMap%Shapes(st2)%bigL, LMAX, &
871 >            plm_j, dlm_j)
872 >
873 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
874 >            CHEBYSHEV_TN, tm_j, dtm_j)
875 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
876              CHEBYSHEV_UN, um_j, dum_j)
877 <      
877 >
878         sigma_j = 0.0d0
879         s_j = 0.0d0
880         eps_j = 0.0d0
# Line 876 | Line 920 | contains  
920               dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
921               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
922            endif
879 –
880 –          sigma_j = sigma_j + plm_j(l,m)*Phunc
881 –          
882 –          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
883 –               Phunc * dlm_j(l,m) * dctjdx
884 –          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
885 –               Phunc * dlm_j(l,m) * dctjdy
886 –          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
887 –               Phunc * dlm_j(l,m) * dctjdz
888 –          
889 –          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
890 –               Phunc * dlm_j(l,m) * dctjdux
891 –          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
892 –               Phunc * dlm_j(l,m) * dctjduy
893 –          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
894 –               Phunc * dlm_j(l,m) * dctjduz
923  
924 +          sigma_j = sigma_j + plm_j(m,l)*Phunc
925 +
926 +          dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
927 +               Phunc * dlm_j(m,l) * dctjdx
928 +          dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
929 +               Phunc * dlm_j(m,l) * dctjdy
930 +          dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
931 +               Phunc * dlm_j(m,l) * dctjdz
932 +
933 +          dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
934 +               Phunc * dlm_j(m,l) * dctjdux
935 +          dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
936 +               Phunc * dlm_j(m,l) * dctjduy
937 +          dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
938 +               Phunc * dlm_j(m,l) * dctjduz
939 +
940         end do
941  
942         do lm = 1, ShapeMap%Shapes(st2)%nRangeFuncs
# Line 919 | Line 963 | contains  
963               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
964            endif
965  
966 <          s_j = s_j + plm_j(l,m)*Phunc
923 <          
924 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
925 <               Phunc * dlm_j(l,m) * dctjdx
926 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
927 <               Phunc * dlm_j(l,m) * dctjdy
928 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
929 <               Phunc * dlm_j(l,m) * dctjdz
930 <          
931 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
932 <               Phunc * dlm_j(l,m) * dctjdux
933 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
934 <               Phunc * dlm_j(l,m) * dctjduy
935 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
936 <               Phunc * dlm_j(l,m) * dctjduz
966 >          s_j = s_j + plm_j(m,l)*Phunc
967  
968 +          dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
969 +               Phunc * dlm_j(m,l) * dctjdx
970 +          dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
971 +               Phunc * dlm_j(m,l) * dctjdy
972 +          dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
973 +               Phunc * dlm_j(m,l) * dctjdz
974 +
975 +          dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
976 +               Phunc * dlm_j(m,l) * dctjdux
977 +          dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
978 +               Phunc * dlm_j(m,l) * dctjduy
979 +          dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
980 +               Phunc * dlm_j(m,l) * dctjduz
981 +
982         end do
983  
984         do lm = 1, ShapeMap%Shapes(st2)%nStrengthFuncs
# Line 961 | Line 1005 | contains  
1005               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1006            endif
1007  
1008 <          eps_j = eps_j + plm_j(l,m)*Phunc
965 <          
966 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
967 <               Phunc * dlm_j(l,m) * dctjdx
968 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
969 <               Phunc * dlm_j(l,m) * dctjdy
970 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
971 <               Phunc * dlm_j(l,m) * dctjdz
972 <          
973 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
974 <               Phunc * dlm_j(l,m) * dctjdux
975 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
976 <               Phunc * dlm_j(l,m) * dctjduy
977 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
978 <               Phunc * dlm_j(l,m) * dctjduz
1008 >          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1009  
1010 +          eps_j = eps_j + plm_j(m,l)*Phunc
1011 +
1012 +          depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1013 +               Phunc * dlm_j(m,l) * dctjdx
1014 +          depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1015 +               Phunc * dlm_j(m,l) * dctjdy
1016 +          depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1017 +               Phunc * dlm_j(m,l) * dctjdz
1018 +
1019 +          depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1020 +               Phunc * dlm_j(m,l) * dctjdux
1021 +          depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
1022 +               Phunc * dlm_j(m,l) * dctjduy
1023 +          depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
1024 +               Phunc * dlm_j(m,l) * dctjduz
1025 +
1026         end do
1027  
1028      endif
# Line 1030 | Line 1076 | contains  
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 <    
1079 >
1080 > !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1081 > !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1082 > !!$
1083 > !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1084 > !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1085 > !!$
1086 > !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1087 >
1088      rtdenom = rij-sigma+s
1089      rt = s / rtdenom
1090  
# Line 1046 | Line 1100 | contains  
1100      drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1101      drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1102      drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1103 <    
1103 >
1104      rt2 = rt*rt
1105      rt3 = rt2*rt
1106      rt5 = rt2*rt3
# Line 1055 | Line 1109 | contains  
1109      rt12 = rt6*rt6
1110      rt126 = rt12 - rt6
1111  
1112 +    pot_temp = 4.0d0 * eps * rt126
1113 +
1114 +    vpair = vpair + pot_temp
1115      if (do_pot) then
1116   #ifdef IS_MPI
1117 <       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1118 <       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1117 >       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1118 >       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1119   #else
1120 <       pot = pot + 4.0d0*eps*rt126*sw
1120 >       pot = pot + pot_temp*sw
1121   #endif
1122      endif
1123 <    
1123 >
1124 > !!$    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
# Line 1081 | Line 1140 | contains  
1140      ! do the torques first since they are easy:
1141      ! remember that these are still in the body fixed axes
1142  
1084 –    txi = dvduxi * sw
1085 –    tyi = dvduyi * sw
1086 –    tzi = dvduzi * sw
1143  
1144 <    txj = dvduxj * sw
1145 <    tyj = dvduyj * sw
1146 <    tzj = dvduzj * sw
1144 > !!$    write(*,*) 'sw = ', sw
1145 > !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1146 > !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1147 > !!$
1148 >    txi =  (dvduzi - dvduyi) * sw
1149 >    tyi =  (dvduxi - dvduzi) * sw
1150 >    tzi =  (dvduyi - dvduxi) * sw
1151  
1152 +    txj = (dvduzj - dvduyj) * sw
1153 +    tyj = (dvduxj - dvduzj) * sw
1154 +    tzj = (dvduyj - dvduxj) * sw
1155 +
1156 + !!$    txi = -dvduxi * sw
1157 + !!$    tyi = -dvduyi * sw
1158 + !!$    tzi = -dvduzi * sw
1159 + !!$
1160 + !!$    txj = dvduxj * sw
1161 + !!$    tyj = dvduyj * sw
1162 + !!$    tzj = dvduzj * sw
1163 +
1164 +    write(*,*) 't1 = ', txi, tyi, tzi
1165 +    write(*,*) 't2 = ', txj, tyj, tzj
1166 +
1167      ! go back to lab frame using transpose of rotation matrix:
1168 <    
1168 >
1169   #ifdef IS_MPI
1170      t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
1171           a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
# Line 1098 | Line 1173 | contains  
1173           a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
1174      t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
1175           a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
1176 <    
1176 >
1177      t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
1178           a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
1179      t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
1180 <            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1180 >         a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1181      t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
1182           a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
1183   #else
1184      t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
1185      t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
1186      t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
1187 <    
1187 >
1188      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1189      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1190      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1191   #endif
1192      ! Now, on to the forces:
1193 <    
1193 >
1194      ! first rotate the i terms back into the lab frame:
1195 <    
1195 >
1196      fxi = dvdxi * sw
1197      fyi = dvdyi * sw
1198      fzi = dvdzi * sw
# Line 1138 | Line 1213 | contains  
1213      fxii = a(1,atom1)*fxi + a(4,atom1)*fyi + a(7,atom1)*fzi
1214      fyii = a(2,atom1)*fxi + a(5,atom1)*fyi + a(8,atom1)*fzi
1215      fzii = a(3,atom1)*fxi + a(6,atom1)*fyi + a(9,atom1)*fzi
1216 <    
1216 >
1217      fxjj = a(1,atom2)*fxj + a(4,atom2)*fyj + a(7,atom2)*fzj
1218      fyjj = a(2,atom2)*fxj + a(5,atom2)*fyj + a(8,atom2)*fzj
1219      fzjj = a(3,atom2)*fxj + a(6,atom2)*fyj + a(9,atom2)*fzj
# Line 1147 | Line 1222 | contains  
1222      fxij = -fxii
1223      fyij = -fyii
1224      fzij = -fzii
1225 <    
1225 >
1226      fxji = -fxjj
1227      fyji = -fyjj
1228      fzji = -fzjj
1229  
1230 <    fxradial = fxii + fxji
1231 <    fyradial = fyii + fyji
1232 <    fzradial = fzii + fzji
1230 >    fxradial = 0.5_dp * (fxii + fxji)
1231 >    fyradial = 0.5_dp * (fyii + fyji)
1232 >    fzradial = 0.5_dp * (fzii + fzji)
1233  
1234   #ifdef IS_MPI
1235      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1236      f_Row(2,atom1) = f_Row(2,atom1) + fyradial
1237      f_Row(3,atom1) = f_Row(3,atom1) + fzradial
1238 <    
1238 >
1239      f_Col(1,atom2) = f_Col(1,atom2) - fxradial
1240      f_Col(2,atom2) = f_Col(2,atom2) - fyradial
1241      f_Col(3,atom2) = f_Col(3,atom2) - fzradial
# Line 1168 | Line 1243 | contains  
1243      f(1,atom1) = f(1,atom1) + fxradial
1244      f(2,atom1) = f(2,atom1) + fyradial
1245      f(3,atom1) = f(3,atom1) + fzradial
1246 <    
1246 >
1247      f(1,atom2) = f(1,atom2) - fxradial
1248      f(2,atom2) = f(2,atom2) - fyradial
1249      f(3,atom2) = f(3,atom2) - fzradial
# Line 1181 | Line 1256 | contains  
1256      id1 = atom1
1257      id2 = atom2
1258   #endif
1259 <    
1259 >
1260      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1261 <      
1261 >
1262         fpair(1) = fpair(1) + fxradial
1263         fpair(2) = fpair(2) + fyradial
1264         fpair(3) = fpair(3) + fzradial
1265 <      
1265 >
1266      endif
1267 <    
1267 >
1268    end subroutine do_shape_pair
1269 <    
1270 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1271 <    
1269 >
1270 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1271 >
1272      ! Purpose: Compute the associated Legendre functions
1273      !          Plm(x) and their derivatives Plm'(x)
1274      ! Input :  x  --- Argument of Plm(x)
# Line 1209 | Line 1284 | contains  
1284      !
1285      ! The original Fortran77 codes can be found here:
1286      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1287 <    
1288 <    real (kind=8), intent(in) :: x
1287 >
1288 >    real (kind=dp), intent(in) :: x
1289      integer, intent(in) :: l, m, lmax
1290 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1290 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1291      integer :: i, j, ls
1292 <    real (kind=8) :: xq, xs
1292 >    real (kind=dp) :: xq, xs
1293  
1294      ! zero out both arrays:
1295      DO I = 0, m
1296         DO J = 0, l
1297 <          PLM(J,I) = 0.0D0
1298 <          DLM(J,I) = 0.0D0
1297 >          PLM(J,I) = 0.0_dp
1298 >          DLM(J,I) = 0.0_dp
1299         end DO
1300      end DO
1301  
1302      ! start with 0,0:
1303      PLM(0,0) = 1.0D0
1304 <  
1304 >
1305      ! x = +/- 1 functions are easy:
1306      IF (abs(X).EQ.1.0D0) THEN
1307         DO I = 1, m
# Line 1253 | Line 1328 | contains  
1328      DO I = 1, l
1329         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1330      enddo
1331 <    
1331 >
1332      DO I = 0, l
1333         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1334      enddo
1335 <    
1335 >
1336      DO I = 0, l
1337         DO J = I+2, m
1338            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1339                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1340         end DO
1341      end DO
1342 <  
1342 >
1343      DLM(0, 0)=0.0D0
1269 –    
1344      DO J = 1, m
1345         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1346      end DO
1347 <    
1347 >
1348      DO I = 1, l
1349         DO J = I, m
1350            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1351         end DO
1352      end DO
1353 <    
1353 >
1354      RETURN
1355    END SUBROUTINE Associated_Legendre
1356  
1357  
1358 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1359 <  
1358 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1359 >
1360      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1361      !          or Ln(x) or Hn(x), and their derivatives
1362      ! Input :  function_type --- Function code
# Line 1301 | Line 1375 | contains  
1375      !
1376      ! The original Fortran77 codes can be found here:
1377      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1378 <  
1378 >
1379      real(kind=8), intent(in) :: x
1380 <    integer, intent(in):: m
1380 >    integer, intent(in):: m, mmax
1381      integer, intent(in):: function_type
1382 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1383 <  
1382 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1383 >
1384      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1385      integer :: k
1386  
# Line 1349 | Line 1423 | contains  
1423         DY0 = DY1
1424         DY1 = DYN
1425      end DO
1426 +
1427 +
1428      RETURN
1429 <    
1429 >
1430    end subroutine Orthogonal_Polynomial
1355 –  
1356 – end module shapes
1431  
1432 < subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1433 <     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1360 <     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1361 <     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1362 <     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1363 <     myAtid, status)
1432 >  subroutine deallocateShapes(this)
1433 >    type(Shape), pointer :: this
1434  
1435 <  use definitions
1436 <  use shapes, only: newShapeType
1437 <  
1438 <  integer :: nContactFuncs
1369 <  integer :: nRangeFuncs
1370 <  integer :: nStrengthFuncs
1371 <  integer :: status
1372 <  integer :: myAtid
1373 <  
1374 <  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1375 <  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1376 <  integer, dimension(nContactFuncs) :: ContactFunctionType        
1377 <  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1378 <  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1379 <  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1380 <  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1381 <  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1382 <  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1383 <  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1384 <  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1385 <  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1386 <  
1387 <  call newShapeType(nContactFuncs, ContactFuncLValue, &
1388 <       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1389 <       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1390 <       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1391 <       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1392 <       myAtid, status)
1435 >    if (associated( this%ContactFuncLValue)) then
1436 >       deallocate(this%ContactFuncLValue)
1437 >       this%ContactFuncLValue => null()
1438 >    end if
1439  
1440 <  return
1441 < end subroutine makeShape
1440 >    if (associated( this%ContactFuncMValue)) then
1441 >       deallocate( this%ContactFuncMValue)
1442 >       this%ContactFuncMValue => null()
1443 >    end if
1444 >    if (associated( this%ContactFunctionType)) then
1445 >       deallocate(this%ContactFunctionType)
1446 >       this%ContactFunctionType => null()
1447 >    end if
1448 >
1449 >    if (associated( this%ContactFuncCoefficient)) then
1450 >       deallocate(this%ContactFuncCoefficient)
1451 >       this%ContactFuncCoefficient => null()
1452 >    end if
1453 >
1454 >    if (associated( this%RangeFuncLValue)) then
1455 >       deallocate(this%RangeFuncLValue)
1456 >       this%RangeFuncLValue => null()
1457 >    end if
1458 >    if (associated( this%RangeFuncMValue)) then
1459 >       deallocate( this%RangeFuncMValue)
1460 >       this%RangeFuncMValue => null()
1461 >    end if
1462 >
1463 >    if (associated( this%RangeFunctionType)) then
1464 >       deallocate( this%RangeFunctionType)
1465 >       this%RangeFunctionType => null()
1466 >    end if
1467 >    if (associated( this%RangeFuncCoefficient)) then
1468 >       deallocate(this%RangeFuncCoefficient)
1469 >       this%RangeFuncCoefficient => null()
1470 >    end if
1471 >
1472 >    if (associated( this%StrengthFuncLValue)) then
1473 >       deallocate(this%StrengthFuncLValue)
1474 >       this%StrengthFuncLValue => null()
1475 >    end if
1476 >
1477 >    if (associated( this%StrengthFuncMValue )) then
1478 >       deallocate(this%StrengthFuncMValue)
1479 >       this%StrengthFuncMValue => null()
1480 >    end if
1481 >
1482 >    if(associated( this%StrengthFunctionType)) then
1483 >       deallocate(this%StrengthFunctionType)
1484 >       this%StrengthFunctionType => null()
1485 >    end if
1486 >    if (associated( this%StrengthFuncCoefficient )) then
1487 >       deallocate(this%StrengthFuncCoefficient)
1488 >       this%StrengthFuncCoefficient => null()
1489 >    end if
1490 >  end subroutine deallocateShapes
1491 >
1492 >  subroutine destroyShapeTypes
1493 >    integer :: i
1494 >    type(Shape), pointer :: thisShape
1495 >
1496 >    ! First walk through and kill the shape
1497 >    do i = 1,ShapeMap%n_shapes
1498 >       thisShape => ShapeMap%Shapes(i)
1499 >       call deallocateShapes(thisShape)
1500 >    end do
1501 >
1502 >    ! set shape map to starting values
1503 >    ShapeMap%n_shapes = 0
1504 >    ShapeMap%currentShape = 0
1505 >
1506 >    if (associated(ShapeMap%Shapes)) then
1507 >       deallocate(ShapeMap%Shapes)
1508 >       ShapeMap%Shapes => null()
1509 >    end if
1510 >
1511 >    if (associated(ShapeMap%atidToShape)) then
1512 >       deallocate(ShapeMap%atidToShape)
1513 >       ShapeMap%atidToShape => null()
1514 >    end if
1515 >
1516 >
1517 >  end subroutine destroyShapeTypes
1518 >
1519 >
1520 > end module shapes

Diff Legend

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