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 1689 by chrisfen, Fri Oct 29 22:28:45 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 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()
102    end type ShapeList
103 <  
103 >
104    type(ShapeList), save :: ShapeMap
105  
106    integer :: lmax
65 –  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
107  
108   contains  
109  
# Line 105 | 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 <      
156 >
157         allocate(ShapeMap%atidToShape(0:ntypes))
158      end if
159 <    
159 >
160      ShapeMap%currentShape = ShapeMap%currentShape + 1
161      current = ShapeMap%currentShape
162  
# Line 150 | 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 188 | 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)
# Line 292 | Line 331 | contains  
331      return
332  
333    end subroutine allocateShape
334 <    
334 >
335    subroutine complete_Shape_FF(status)
336      integer :: status
337      integer :: i, j, l, m, lm, function_type
# Line 306 | Line 345 | contains  
345         status = -1
346         return
347      end if
348 <    
348 >
349      nAtypes = getSize(atypes)
350  
351      if (nAtypes == 0) then
# Line 316 | Line 355 | contains  
355  
356      ! atypes comes from c side
357      do i = 0, nAtypes
358 <      
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 332 | Line 371 | contains  
371  
372            ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
373            ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
374 <          
374 >
375         endif
376 <      
376 >
377      end do
378  
379      haveShapeMap = .true.
380 <    
380 >
381    end subroutine complete_Shape_FF
382 <    
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  
# Line 351 | Line 390 | contains  
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 362 | 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 380 | 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 411 | 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 438 | 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(LMAX,MMAX), dlm_i(LMAX,MMAX)
484 <    real (kind=dp) :: plm_j(LMAX,MMAX), dlm_j(LMAX,MMAX)
485 <    real (kind=dp) :: tm_i(MMAX), dtm_i(MMAX), um_i(MMAX), dum_i(MMAX)
486 <    real (kind=dp) :: tm_j(MMAX), dtm_j(MMAX), um_j(MMAX), dum_j(MMAX)
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  
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 461 | 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 503 | 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  
563         xi2 = xi*xi
564         yi2 = yi*yi
565 <       zi2 = zi*zi
524 <      
525 <       proji = sqrt(xi2 + yi2)
526 <       proji3 = proji*proji*proji
527 <      
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
538 –       dcpidx = 1.0d0 / proji - xi2 / proji3
539 –       dcpidy = - xi * yi / proji3
606         dcpidz = 0.0d0
607 <       dcpidux = xi * yi * zi / proji3
608 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
543 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
544 <      
607 >       dcpiduz = 0.0d0
608 >
609         spi = yi / proji
546 –       dspidx = - xi * yi / proji3
547 –       dspidy = 1.0d0 / proji - yi2 / proji3
610         dspidz = 0.0d0
611 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
550 <       dspiduy = xi * yi * zi / proji3
551 <       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, &
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, MMAX, &
618              CHEBYSHEV_TN, tm_i, dtm_i)
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 588 | Line 648 | contains  
648            function_type = ShapeMap%Shapes(st1)%ContactFunctionType(lm)
649  
650            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
591 –           !  write(*,*) tm_i(m), ' is tm_i'
651               Phunc = coeff * tm_i(m)
652               dPhuncdX = coeff * dtm_i(m) * dcpidx
653               dPhuncdY = coeff * dtm_i(m) * dcpidy
# Line 606 | 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
610 <          write(*,*) plm_i(l,m), l, m
611 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
612 <               Phunc * dlm_i(l,m) * dctidx
613 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
614 <               Phunc * dlm_i(l,m) * dctidy
615 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
616 <               Phunc * dlm_i(l,m) * dctidz
617 <          
618 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
619 <               Phunc * dlm_i(l,m) * dctidux
620 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
621 <               Phunc * dlm_i(l,m) * dctiduy
622 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
623 <               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 629 | 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 648 | 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
652 <          
653 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
654 <               Phunc * dlm_i(l,m) * dctidx
655 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
656 <               Phunc * dlm_i(l,m) * dctidy
657 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
658 <               Phunc * dlm_i(l,m) * dctidz
659 <          
660 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
661 <               Phunc * dlm_i(l,m) * dctidux
662 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
663 <               Phunc * dlm_i(l,m) * dctiduy
664 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
665 <               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 689 | Line 748 | contains  
748               dPhuncdUy = coeff*(spi * dum_i(m-1)*dcpiduy + dspiduy *um_i(m-1))
749               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
750            endif
692 –          !write(*,*) eps_i, plm_i(l,m), Phunc
693 –          eps_i = eps_i + plm_i(l,m)*Phunc
694 –          
695 –          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
696 –               Phunc * dlm_i(l,m) * dctidx
697 –          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
698 –               Phunc * dlm_i(l,m) * dctidy
699 –          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
700 –               Phunc * dlm_i(l,m) * dctidz
701 –          
702 –          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
703 –               Phunc * dlm_i(l,m) * dctidux
704 –          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
705 –               Phunc * dlm_i(l,m) * dctiduy
706 –          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
707 –               Phunc * dlm_i(l,m) * dctiduz      
751  
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
712 –      
713 –       ! 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 735 | 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 748 | 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
760 –      
761 –       projj = sqrt(xj2 + yj2)
762 –       projj3 = projj*projj*projj
763 –      
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 <      
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
773 –       dcpjdx = 1.0d0 / projj - xj2 / projj3
774 –       dcpjdy = - xj * yj / projj3
859         dcpjdz = 0.0d0
860 <       dcpjdux = xj * yj * zj / projj3
861 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
778 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
779 <      
860 >       dcpjduz = 0.0d0
861 >
862         spj = yj / projj
781 –       dspjdx = - xj * yj / projj3
782 –       dspjdy = 1.0d0 / projj - yj2 / projj3
863         dspjdz = 0.0d0
864 <       dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
865 <       dspjduy = xj * yj * zj / projj3
866 <       dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
867 <      
868 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
869 <            ShapeMap%Shapes(st2)%bigM, LMAX, &
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 <      
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 839 | 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
842 –
843 –          sigma_j = sigma_j + plm_j(l,m)*Phunc
844 –          
845 –          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
846 –               Phunc * dlm_j(l,m) * dctjdx
847 –          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
848 –               Phunc * dlm_j(l,m) * dctjdy
849 –          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
850 –               Phunc * dlm_j(l,m) * dctjdz
851 –          
852 –          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
853 –               Phunc * dlm_j(l,m) * dctjdux
854 –          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
855 –               Phunc * dlm_j(l,m) * dctjduy
856 –          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
857 –               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 882 | 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
886 <          
887 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
888 <               Phunc * dlm_j(l,m) * dctjdx
889 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
890 <               Phunc * dlm_j(l,m) * dctjdy
891 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
892 <               Phunc * dlm_j(l,m) * dctjdz
893 <          
894 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
895 <               Phunc * dlm_j(l,m) * dctjdux
896 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
897 <               Phunc * dlm_j(l,m) * dctjduy
898 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
899 <               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 924 | 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
928 <          
929 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
930 <               Phunc * dlm_j(l,m) * dctjdx
931 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
932 <               Phunc * dlm_j(l,m) * dctjdy
933 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
934 <               Phunc * dlm_j(l,m) * dctjdz
935 <          
936 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
937 <               Phunc * dlm_j(l,m) * dctjdux
938 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
939 <               Phunc * dlm_j(l,m) * dctjduy
940 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
941 <               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 977 | Line 1060 | contains  
1060      dsduxj = 0.5*dsjdux
1061      dsduyj = 0.5*dsjduy
1062      dsduzj = 0.5*dsjduz
1063 <    !write(*,*) eps_i, eps_j
1063 >
1064      eps = sqrt(eps_i * eps_j)
1065  
1066      depsdxi = eps_j * depsidx / (2.0d0 * eps)
# Line 993 | 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 1009 | 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 1018 | 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 1044 | Line 1140 | contains  
1140      ! do the torques first since they are easy:
1141      ! remember that these are still in the body fixed axes
1142  
1047 –    txi = dvduxi * sw
1048 –    tyi = dvduyi * sw
1049 –    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 1061 | 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 1101 | 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 1110 | 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 1131 | 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 1144 | 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 <    
1269 >
1270    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1271  
1272      ! Purpose: Compute the associated Legendre functions
# Line 1172 | 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 <    
1287 >
1288      real (kind=dp), intent(in) :: x
1289      integer, intent(in) :: l, m, lmax
1290      real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
# Line 1189 | Line 1301 | contains  
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 1244 | Line 1356 | contains  
1356  
1357  
1358    subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1359 <  
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 1263 | 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, mmax
1381      integer, intent(in):: function_type
1382      real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1383 <  
1383 >
1384      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1385      integer :: k
1386  
# Line 1311 | Line 1423 | contains  
1423         DY0 = DY1
1424         DY1 = DYN
1425      end DO
1426 +
1427 +
1428      RETURN
1429 <    
1429 >
1430    end subroutine Orthogonal_Polynomial
1317 –  
1318 – end module shapes
1431  
1432 < subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1433 <     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1322 <     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1323 <     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1324 <     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1325 <     myATID, status)
1432 >  subroutine deallocateShapes(this)
1433 >    type(Shape), pointer :: this
1434  
1435 <  use definitions
1436 <  use shapes, only: newShapeType
1437 <  
1438 <  integer :: nContactFuncs
1331 <  integer :: nRangeFuncs
1332 <  integer :: nStrengthFuncs
1333 <  integer :: status
1334 <  integer :: myATID
1335 <  
1336 <  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1337 <  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1338 <  integer, dimension(nContactFuncs) :: ContactFunctionType        
1339 <  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1340 <  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1341 <  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1342 <  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1343 <  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1344 <  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1345 <  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1346 <  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1347 <  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1348 <  
1349 <  call newShapeType(nContactFuncs, ContactFuncLValue, &
1350 <       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1351 <       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1352 <       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1353 <       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1354 <       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 < subroutine completeShapeFF(status)
1449 >    if (associated( this%ContactFuncCoefficient)) then
1450 >       deallocate(this%ContactFuncCoefficient)
1451 >       this%ContactFuncCoefficient => null()
1452 >    end if
1453  
1454 <  use shapes, only: complete_Shape_FF
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 <  integer, intent(out)  :: status
1464 <  integer :: myStatus
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 <  myStatus = 0
1472 >    if (associated( this%StrengthFuncLValue)) then
1473 >       deallocate(this%StrengthFuncLValue)
1474 >       this%StrengthFuncLValue => null()
1475 >    end if
1476  
1477 <  call complete_Shape_FF(myStatus)
1477 >    if (associated( this%StrengthFuncMValue )) then
1478 >       deallocate(this%StrengthFuncMValue)
1479 >       this%StrengthFuncMValue => null()
1480 >    end if
1481  
1482 <  status = myStatus
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 <  return
1493 < end subroutine completeShapeFF
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