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 1647 by chrisfen, Tue Oct 26 18:03:12 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 <  real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
66 <  real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
67 <  real (kind=dp), allocatable, dimension(:) :: tm_j, dtm_j, um_j, dum_j
68 <
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 :: myAtid
122 >    integer :: c_ident
123 >    integer :: myATID
124      integer :: bigL
125      integer :: bigM
126      integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
# Line 105 | 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 <      
157 >
158         allocate(ShapeMap%atidToShape(ntypes))
159      end if
160 <    
160 >
161      ShapeMap%currentShape = ShapeMap%currentShape + 1
162      current = ShapeMap%currentShape
163  
# Line 128 | Line 168 | contains  
168         return
169      endif
170  
171 <    call getElementProperty(atypes, myAtid, "c_ident", me)
172 <    ShapeMap%atidToShape(me)                         = current
173 <    ShapeMap%Shapes(current)%atid                    = me
171 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
172 >
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 149 | 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 187 | 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)
235      endif
# Line 253 | Line 295 | contains  
295         stat = -1
296         return
297      endif
298 <    
298 >
299      if (associated(myShape%strengthFuncLValue)) then
300         deallocate(myShape%strengthFuncLValue)
301      endif
# Line 287 | Line 329 | contains  
329         return
330      endif
331  
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 302 | 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 = 1, nAtypes
359 <      
360 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
361 <      
359 >    
360 >       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
361 >       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
362 >        
363         if (thisProperty) then
318 –          
364            ShapeMap%currentShape = ShapeMap%currentShape + 1
365            current = ShapeMap%currentShape
366  
367 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
368 <          ShapeMap%atidToShape(thisIP) = current
324 <          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 >
393      integer, intent(in) :: atom1, atom2
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 354 | 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 372 | 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 403 | 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 430 | 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)
492 +    real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
493 +
494      if (.not.haveShapeMap) then
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.
440 –        
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 459 | Line 525 | contains  
525   #endif
526  
527      ! use the atid to find the shape type (st) for each atom:
462 –
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 +
535         sigma_i = ShapeMap%Shapes(st1)%sigma
536         s_i = ShapeMap%Shapes(st1)%sigma
537         eps_i = ShapeMap%Shapes(st1)%epsilon
# Line 490 | 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)
573 <      
573 >
574   #endif
575 <      
575 >       xihat = xi / rij
576 >       yihat = yi / rij
577 >       zihat = zi / rij
578         xi2 = xi*xi
579         yi2 = yi*yi
580 <       zi2 = zi*zi
511 <      
512 <       proji = sqrt(xi2 + yi2)
513 <       proji3 = proji*proji*proji
514 <      
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
524 –       dcpidx = 1.0d0 / proji - xi2 / proji3
525 –       dcpidy = - xi * yi / proji3
621         dcpidz = 0.0d0
622 <       dcpidux = xi * yi * zi / proji3
623 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
529 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
530 <      
622 >       dcpiduz = 0.0d0
623 >
624         spi = yi / proji
532 –       dspidx = - xi * yi / proji3
533 –       dspidy = 1.0d0 / proji - yi2 / proji3
625         dspidz = 0.0d0
626 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
536 <       dspiduy = xi * yi * zi / proji3
537 <       dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
626 >       dspiduz = 0.0d0
627  
628 <       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
629 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
628 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
629 >            ShapeMap%Shapes(st1)%bigL, LMAX, &
630 >            plm_i, dlm_i)
631  
632 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
632 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
633              CHEBYSHEV_TN, tm_i, dtm_i)
634 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
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 577 | 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 590 | Line 680 | contains  
680               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
681            endif
682  
683 <          sigma_i = sigma_i + plm_i(l,m)*Phunc
684 <          
685 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
686 <               Phunc * dlm_i(l,m) * dctidx
687 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
688 <               Phunc * dlm_i(l,m) * dctidy
689 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
690 <               Phunc * dlm_i(l,m) * dctidz
691 <          
692 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
693 <               Phunc * dlm_i(l,m) * dctidux
694 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
695 <               Phunc * dlm_i(l,m) * dctiduy
696 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
697 <               Phunc * dlm_i(l,m) * dctiduz
698 <
683 >          sigma_i = sigma_i + plm_i(m,l)*Phunc
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
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 > !!$          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 613 | 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 632 | Line 724 | contains  
724               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
725            endif
726  
727 <          s_i = s_i + plm_i(l,m)*Phunc
636 <          
637 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
638 <               Phunc * dlm_i(l,m) * dctidx
639 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
640 <               Phunc * dlm_i(l,m) * dctidy
641 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
642 <               Phunc * dlm_i(l,m) * dctidz
643 <          
644 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
645 <               Phunc * dlm_i(l,m) * dctidux
646 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
647 <               Phunc * dlm_i(l,m) * dctiduy
648 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
649 <               Phunc * dlm_i(l,m) * dctiduz      
727 >          s_i = s_i + plm_i(m,l)*Phunc
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 +
736 +          dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
737 +               Phunc * dlm_i(m,l) * dctidux
738 +          dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
739 +               Phunc * dlm_i(m,l) * dctiduy
740 +          dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
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 674 | Line 766 | contains  
766               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
767            endif
768  
769 <          eps_i = eps_i + plm_i(l,m)*Phunc
770 <          
771 <          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
772 <               Phunc * dlm_i(l,m) * dctidx
773 <          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
774 <               Phunc * dlm_i(l,m) * dctidy
775 <          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
776 <               Phunc * dlm_i(l,m) * dctidz
777 <          
778 <          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
779 <               Phunc * dlm_i(l,m) * dctidux
780 <          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
781 <               Phunc * dlm_i(l,m) * dctiduy
782 <          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
783 <               Phunc * dlm_i(l,m) * dctiduz      
769 >          eps_i = eps_i + plm_i(m,l)*Phunc
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 >
778 >          depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
779 >               Phunc * dlm_i(m,l) * dctidux
780 >          depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
781 >               Phunc * dlm_i(m,l) * dctiduy
782 >          depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
783 >               Phunc * dlm_i(m,l) * dctiduz      
784  
785         end do
786  
787      endif
696 –      
697 –       ! 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 719 | 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 732 | 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
744 –      
745 –       projj = sqrt(xj2 + yj2)
746 –       projj3 = projj*projj*projj
747 –      
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
757 –       dcpjdx = 1.0d0 / projj - xj2 / projj3
758 –       dcpjdy = - xj * yj / projj3
879         dcpjdz = 0.0d0
880 <       dcpjdux = xj * yj * zj / projj3
881 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
762 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
763 <      
880 >       dcpjduz = 0.0d0
881 >
882         spj = yj / projj
765 –       dspjdx = - xj * yj / projj3
766 –       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 <      
888 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
889 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
890 <      
891 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
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 >
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, &
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 810 | 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 822 | 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
825 –
826 –          sigma_j = sigma_j + plm_j(l,m)*Phunc
827 –          
828 –          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
829 –               Phunc * dlm_j(l,m) * dctjdx
830 –          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
831 –               Phunc * dlm_j(l,m) * dctjdy
832 –          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
833 –               Phunc * dlm_j(l,m) * dctjdz
834 –          
835 –          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
836 –               Phunc * dlm_j(l,m) * dctjdux
837 –          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
838 –               Phunc * dlm_j(l,m) * dctjduy
839 –          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
840 –               Phunc * dlm_j(l,m) * dctjduz
943  
944 +          sigma_j = sigma_j + plm_j(m,l)*Phunc
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 +
953 +          dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
954 +               Phunc * dlm_j(m,l) * dctjdux
955 +          dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
956 +               Phunc * dlm_j(m,l) * dctjduy
957 +          dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
958 +               Phunc * dlm_j(m,l) * dctjduz
959 +
960         end do
961  
962         do lm = 1, ShapeMap%Shapes(st2)%nRangeFuncs
# Line 852 | 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 865 | Line 983 | contains  
983               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
984            endif
985  
986 <          s_j = s_j + plm_j(l,m)*Phunc
869 <          
870 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
871 <               Phunc * dlm_j(l,m) * dctjdx
872 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
873 <               Phunc * dlm_j(l,m) * dctjdy
874 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
875 <               Phunc * dlm_j(l,m) * dctjdz
876 <          
877 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
878 <               Phunc * dlm_j(l,m) * dctjdux
879 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
880 <               Phunc * dlm_j(l,m) * dctjduy
881 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
882 <               Phunc * dlm_j(l,m) * dctjduz
986 >          s_j = s_j + plm_j(m,l)*Phunc
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 +
995 +          dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
996 +               Phunc * dlm_j(m,l) * dctjdux
997 +          dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
998 +               Phunc * dlm_j(m,l) * dctjduy
999 +          dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
1000 +               Phunc * dlm_j(m,l) * dctjduz
1001 +
1002         end do
1003  
1004         do lm = 1, ShapeMap%Shapes(st2)%nStrengthFuncs
# Line 907 | Line 1025 | contains  
1025               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1026            endif
1027  
1028 <          eps_j = eps_j + plm_j(l,m)*Phunc
911 <          
912 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
913 <               Phunc * dlm_j(l,m) * dctjdx
914 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
915 <               Phunc * dlm_j(l,m) * dctjdy
916 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
917 <               Phunc * dlm_j(l,m) * dctjdz
918 <          
919 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
920 <               Phunc * dlm_j(l,m) * dctjdux
921 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
922 <               Phunc * dlm_j(l,m) * dctjduy
923 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
924 <               Phunc * dlm_j(l,m) * dctjduz
1028 > !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1029  
1030 +          eps_j = eps_j + plm_j(m,l)*Phunc
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 +
1039 +          depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1040 +               Phunc * dlm_j(m,l) * dctjdux
1041 +          depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
1042 +               Phunc * dlm_j(m,l) * dctjduy
1043 +          depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
1044 +               Phunc * dlm_j(m,l) * dctjduz
1045 +
1046         end do
1047  
1048      endif
# Line 930 | 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 962 | 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 976 | 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 1001 | 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 1023 | 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
1167 >    txi = 0.0d0
1168 >    tyi = 0.0d0
1169 >    tzi = 0.0d0
1170  
1171 <    txj = dvduxj * sw
1172 <    tyj = dvduyj * sw
1173 <    tzj = dvduzj * sw
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 1044 | 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:
1066 –    
1067 –    fxi = dvdxi * sw
1068 –    fyi = dvdyi * sw
1069 –    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 1084 | 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 1093 | 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 1114 | 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 1127 | 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 <    
1299 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1300 <    
1298 >
1299 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1300 >
1301      ! Purpose: Compute the associated Legendre functions
1302      !          Plm(x) and their derivatives Plm'(x)
1303      ! Input :  x  --- Argument of Plm(x)
# Line 1155 | 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 <    
1317 <    real (kind=8), intent(in) :: x
1316 >
1317 >    real (kind=dp), intent(in) :: x
1318      integer, intent(in) :: l, m, lmax
1319 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1319 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1320      integer :: i, j, ls
1321 <    real (kind=8) :: xq, xs
1321 >    real (kind=dp) :: xq, xs
1322  
1323      ! zero out both arrays:
1324      DO I = 0, m
1325         DO J = 0, l
1326 <          PLM(J,I) = 0.0D0
1327 <          DLM(J,I) = 0.0D0
1326 >          PLM(J,I) = 0.0_dp
1327 >          DLM(J,I) = 0.0_dp
1328         end DO
1329      end DO
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 1199 | Line 1357 | contains  
1357      DO I = 1, l
1358         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1359      enddo
1360 <    
1360 >
1361      DO I = 0, l
1362         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1363      enddo
1364 <    
1364 >
1365      DO I = 0, l
1366         DO J = I+2, m
1367            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1368                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1369         end DO
1370      end DO
1371 <  
1371 >
1372      DLM(0, 0)=0.0D0
1215 –    
1373      DO J = 1, m
1374         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1375      end DO
1376 <    
1376 >
1377      DO I = 1, l
1378         DO J = I, m
1379            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1380         end DO
1381      end DO
1382 <    
1382 >
1383      RETURN
1384    END SUBROUTINE Associated_Legendre
1385  
1386  
1387 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1388 <  
1387 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
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 1247 | 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
1409 >    integer, intent(in):: m, mmax
1410      integer, intent(in):: function_type
1411 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1412 <  
1411 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1412 >
1413      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1414      integer :: k
1415  
# Line 1295 | Line 1452 | contains  
1452         DY0 = DY1
1453         DY1 = DYN
1454      end DO
1298 –    RETURN
1299 –    
1300 –  end subroutine Orthogonal_Polynomial
1301 –  
1302 – end module shapes
1455  
1304 – subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1305 –     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1306 –     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1307 –     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1308 –     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1309 –     myAtid, status)
1456  
1457 <  use definitions
1312 <  use shapes, only: newShapeType
1313 <  
1314 <  integer :: nContactFuncs
1315 <  integer :: nRangeFuncs
1316 <  integer :: nStrengthFuncs
1317 <  integer :: status
1318 <  integer :: myAtid
1319 <  
1320 <  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1321 <  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1322 <  integer, dimension(nContactFuncs) :: ContactFunctionType        
1323 <  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1324 <  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1325 <  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1326 <  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1327 <  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1328 <  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1329 <  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1330 <  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1331 <  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1332 <  
1333 <  call newShapeType(nContactFuncs, ContactFuncLValue, &
1334 <       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1335 <       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1336 <       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1337 <       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1338 <       myAtid, status)
1457 >    RETURN
1458  
1459 <  return
1341 < end subroutine makeShape
1459 >  end subroutine Orthogonal_Polynomial
1460  
1461 < subroutine completeShapeFF(status)
1462 <
1345 <  use shapes, only: complete_Shape_FF
1461 >  subroutine deallocateShapes(this)
1462 >    type(Shape), pointer :: this
1463  
1464 <  integer, intent(out)  :: status
1465 <  integer :: myStatus
1464 >    if (associated( this%ContactFuncLValue)) then
1465 >       deallocate(this%ContactFuncLValue)
1466 >       this%ContactFuncLValue => null()
1467 >    end if
1468  
1469 <  myStatus = 0
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 <  call complete_Shape_FF(myStatus)
1478 >    if (associated( this%ContactFuncCoefficient)) then
1479 >       deallocate(this%ContactFuncCoefficient)
1480 >       this%ContactFuncCoefficient => null()
1481 >    end if
1482  
1483 <  status = myStatus
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 <  return
1493 < end subroutine completeShapeFF
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 +    if (associated( this%StrengthFuncLValue)) then
1502 +       deallocate(this%StrengthFuncLValue)
1503 +       this%StrengthFuncLValue => null()
1504 +    end if
1505 +
1506 +    if (associated( this%StrengthFuncMValue )) then
1507 +       deallocate(this%StrengthFuncMValue)
1508 +       this%StrengthFuncMValue => null()
1509 +    end if
1510 +
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 +  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