ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/UseTheForce/DarkSide/shapes.F90
(Generate patch)

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2214 by chrisfen, Wed Apr 27 20:14:03 2005 UTC

# Line 55 | 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 68 | Line 68 | module shapes
68    public :: do_shape_pair
69    public :: newShapeType
70    public :: complete_Shape_FF
71 +  public :: destroyShapeTypes
72  
72 –
73    type, private :: Shape
74       integer :: atid
75       integer :: nContactFuncs
# Line 93 | Line 93 | module shapes
93       real ( kind = dp )  :: epsilon
94       real ( kind = dp )  :: sigma
95    end type Shape
96 <  
96 >
97    type, private :: ShapeList
98       integer :: n_shapes = 0
99       integer :: currentShape = 0
100 <     type (Shape), pointer :: Shapes(:)      => null()
101 <     integer, pointer      :: atidToShape(:) => null()
100 >     type(Shape), pointer :: Shapes(:)      => null()
101 >     integer, pointer     :: atidToShape(:) => null()
102    end type ShapeList
103    
104    type(ShapeList), save :: ShapeMap
105 <
105 >  
106    integer :: lmax
107 <
107 >  
108   contains  
109 <
109 >  
110    subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
111         ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115 <       myATID, status)
116 <
115 >       c_ident, status)
116 >    
117      integer :: nContactFuncs
118      integer :: nRangeFuncs
119      integer :: nStrengthFuncs
120      integer :: shape_ident
121      integer :: status
122 +    integer :: c_ident
123      integer :: myATID
124      integer :: bigL
125      integer :: bigM
# Line 144 | Line 145 | contains  
145  
146         call getMatchingElementList(atypes, "is_Shape", .true., &
147              nShapeTypes, MatchList)
148 <      
148 >
149         call getMatchingElementList(atypes, "is_LennardJones", .true., &
150              nLJTypes, MatchList)
151 <      
151 >
152         ShapeMap%n_shapes = nShapeTypes + nLJTypes
153 <      
153 >
154         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
155 <      
155 >
156         ntypes = getSize(atypes)
157 <      
158 <       allocate(ShapeMap%atidToShape(0:ntypes))
157 >
158 >       allocate(ShapeMap%atidToShape(ntypes))
159      end if
160 <    
160 >
161      ShapeMap%currentShape = ShapeMap%currentShape + 1
162      current = ShapeMap%currentShape
163  
# Line 167 | Line 168 | contains  
168         return
169      endif
170  
171 <    call getElementProperty(atypes, myATID, 'c_ident', me)
171 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
172  
173 <    ShapeMap%atidToShape(me)                         = current
174 <    ShapeMap%Shapes(current)%atid                    = me
173 >    ShapeMap%atidToShape(myATID)                     = current
174 >    ShapeMap%Shapes(current)%atid                    = myATID
175      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
176      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
177      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 189 | 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 227 | Line 228 | contains  
228      type(Shape), intent(inout) :: myShape
229      integer, intent(out) :: stat
230      integer :: alloc_stat
231 <
231 >
232      stat = 0
233      if (associated(myShape%contactFuncLValue)) then
234         deallocate(myShape%contactFuncLValue)
# Line 331 | Line 332 | contains  
332      return
333  
334    end subroutine allocateShape
335 <    
335 >
336    subroutine complete_Shape_FF(status)
337      integer :: status
338      integer :: i, j, l, m, lm, function_type
339      real(kind=dp) :: thisDP, sigma
340 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
340 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
341      logical :: thisProperty
342  
343      status = 0
# Line 345 | Line 346 | contains  
346         status = -1
347         return
348      end if
349 <    
349 >
350      nAtypes = getSize(atypes)
351  
352      if (nAtypes == 0) then
353         status = -1
354         return
355 <    end if
355 >    end if      
356  
357      ! atypes comes from c side
358 <    do i = 0, nAtypes
359 <      
360 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
361 <      
358 >    do i = 1, nAtypes
359 >    
360 >       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
361 >       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
362 >        
363         if (thisProperty) then
362 –          
364            ShapeMap%currentShape = ShapeMap%currentShape + 1
365            current = ShapeMap%currentShape
366  
367 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
368 <          ShapeMap%atidToShape(thisIP) = current
368 <          ShapeMap%Shapes(current)%atid = thisIP
367 >          ShapeMap%atidToShape(myATID) = current
368 >          ShapeMap%Shapes(current)%atid = myATID
369  
370            ShapeMap%Shapes(current)%isLJ = .true.
371  
372 <          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
373 <          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
374 <          
372 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
373 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
374 >
375         endif
376 <      
376 >
377      end do
378  
379      haveShapeMap = .true.
380 <    
380 >
381 > !    do i = 1, ShapeMap%n_shapes
382 > !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
383 > !    end do
384 >
385    end subroutine complete_Shape_FF
386 <    
386 >
387    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
388         pot, A, f, t, do_pot)
389 <    
389 >
390      INTEGER, PARAMETER:: LMAX         = 64
391      INTEGER, PARAMETER:: MMAX         = 64
392  
# Line 453 | 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 480 | Line 484 | contains  
484      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
485      real (kind=dp) :: fxradial, fyradial, fzradial
486  
487 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
488 +
489      real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
490      real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
491      real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
# Line 489 | Line 495 | contains  
495         call handleError("calc_shape", "NO SHAPEMAP!!!!")
496         return      
497      endif
498 <    
498 >
499      !! We assume that the rotation matrices have already been calculated
500      !! and placed in the A array.
495 –
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 516 | Line 527 | contains  
527      ! use the atid to find the shape type (st) for each atom:
528      st1 = ShapeMap%atidToShape(atid1)
529      st2 = ShapeMap%atidToShape(atid2)
530 +    
531 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
532  
533      if (ShapeMap%Shapes(st1)%isLJ) then
534  
# Line 545 | 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)
560 –      
561 – #endif
573  
574 + #endif
575 +       xihat = xi / rij
576 +       yihat = yi / rij
577 +       zihat = zi / rij
578         xi2 = xi*xi
579         yi2 = yi*yi
580         zi2 = zi*zi            
# Line 571 | Line 586 | contains  
586         dctidx = - zi * xi / r3
587         dctidy = - zi * yi / r3
588         dctidz = 1.0d0 / rij - zi2 / r3
589 <       dctidux = - (zi * xi2) / r3
590 <       dctiduy = - (zi * yi2) / r3
591 <       dctiduz = zi / rij - (zi2 * zi) / r3
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:
# Line 601 | Line 616 | contains  
616            dspidux = - (yi * xi2) / proji3
617            dspiduy = yi / proji - (yi2 * yi) / proji3
618         endif
619 <      
619 >
620         cpi = xi / proji
621         dcpidz = 0.0d0
622         dcpiduz = 0.0d0
623 <      
623 >
624         spi = yi / proji
625         dspidz = 0.0d0
626         dspiduz = 0.0d0
# Line 618 | Line 633 | contains  
633              CHEBYSHEV_TN, tm_i, dtm_i)
634         call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
635              CHEBYSHEV_UN, um_i, dum_i)
636 <      
636 >
637         sigma_i = 0.0d0
638         s_i = 0.0d0
639         eps_i = 0.0d0
# Line 652 | 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 666 | Line 681 | contains  
681            endif
682  
683            sigma_i = sigma_i + plm_i(m,l)*Phunc
684 <
684 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux
685 > !!$          write(*,*) 'Phunc = ', Phunc
686            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
687                 Phunc * dlm_i(m,l) * dctidx
688            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
689                 Phunc * dlm_i(m,l) * dctidy
690            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
691                 Phunc * dlm_i(m,l) * dctidz
676 –          
692            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
693                 Phunc * dlm_i(m,l) * dctidux
694            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
695                 Phunc * dlm_i(m,l) * dctiduy
696            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
697                 Phunc * dlm_i(m,l) * dctiduz
698 <
698 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
699 > !!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
700 > !!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
701         end do
702  
703         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 688 | 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 708 | Line 725 | contains  
725            endif
726  
727            s_i = s_i + plm_i(m,l)*Phunc
728 <          
728 >
729            dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
730                 Phunc * dlm_i(m,l) * dctidx
731            dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
732                 Phunc * dlm_i(m,l) * dctidy
733            dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
734                 Phunc * dlm_i(m,l) * dctidz
735 <          
735 >
736            dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
737                 Phunc * dlm_i(m,l) * dctidux
738            dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 724 | Line 741 | contains  
741                 Phunc * dlm_i(m,l) * dctiduz      
742  
743         end do
744 <              
744 >
745         do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
746            l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
747            m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
748            coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
749            function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
750 <          
750 >
751            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
752               Phunc = coeff * tm_i(m)
753               dPhuncdX = coeff * dtm_i(m) * dcpidx
754               dPhuncdY = coeff * dtm_i(m) * dcpidy
755               dPhuncdZ = coeff * dtm_i(m) * dcpidz
756 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
756 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
757               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
758               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
759            else
# Line 750 | Line 767 | contains  
767            endif
768  
769            eps_i = eps_i + plm_i(m,l)*Phunc
770 <          
770 >
771            depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
772                 Phunc * dlm_i(m,l) * dctidx
773            depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
774                 Phunc * dlm_i(m,l) * dctidy
775            depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
776                 Phunc * dlm_i(m,l) * dctidz
777 <          
777 >
778            depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
779                 Phunc * dlm_i(m,l) * dctidux
780            depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 768 | Line 785 | contains  
785         end do
786  
787      endif
771 –      
772 –       ! 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 794 | 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 807 | 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
839         ctj = zj / rij
840 <      
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 = - (zi * xj2) / r3
848 <       dctjduy = - (zj * yj2) / r3
849 <       dctjduz = zj / rij - (zj2 * zj) / r3
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  
# Line 858 | Line 878 | contains  
878         cpj = xj / projj
879         dcpjdz = 0.0d0
880         dcpjduz = 0.0d0
881 <      
881 >
882         spj = yj / projj
883         dspjdz = 0.0d0
884         dspjduz = 0.0d0
885  
886  
887 <       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
888 <       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
887 > !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
888 > !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
889         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
890              ShapeMap%Shapes(st2)%bigL, LMAX, &
891              plm_j, dlm_j)
892 <      
892 >
893         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
894              CHEBYSHEV_TN, tm_j, dtm_j)
895         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
896              CHEBYSHEV_UN, um_j, dum_j)
897 <      
897 >
898         sigma_j = 0.0d0
899         s_j = 0.0d0
900         eps_j = 0.0d0
# Line 908 | 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 920 | Line 940 | contains  
940               dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
941               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
942            endif
943 <
943 >
944            sigma_j = sigma_j + plm_j(m,l)*Phunc
945 <          
945 >
946            dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
947                 Phunc * dlm_j(m,l) * dctjdx
948            dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
949                 Phunc * dlm_j(m,l) * dctjdy
950            dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
951                 Phunc * dlm_j(m,l) * dctjdz
952 <          
952 >
953            dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
954                 Phunc * dlm_j(m,l) * dctjdux
955            dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
# Line 950 | 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 964 | Line 984 | contains  
984            endif
985  
986            s_j = s_j + plm_j(m,l)*Phunc
987 <          
987 >
988            dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
989                 Phunc * dlm_j(m,l) * dctjdx
990            dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
991                 Phunc * dlm_j(m,l) * dctjdy
992            dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
993                 Phunc * dlm_j(m,l) * dctjdz
994 <          
994 >
995            dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
996                 Phunc * dlm_j(m,l) * dctjdux
997            dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 1005 | Line 1025 | contains  
1025               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1026            endif
1027  
1028 <          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1028 > !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1029  
1030            eps_j = eps_j + plm_j(m,l)*Phunc
1031 <          
1031 >
1032            depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1033                 Phunc * dlm_j(m,l) * dctjdx
1034            depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1035                 Phunc * dlm_j(m,l) * dctjdy
1036            depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1037                 Phunc * dlm_j(m,l) * dctjdz
1038 <          
1038 >
1039            depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1040                 Phunc * dlm_j(m,l) * dctjdux
1041            depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 1030 | 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 1062 | 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 1076 | 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
1082 – !!$
1083 – !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1105   !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1106   !!$
1107   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
# Line 1088 | Line 1109 | contains  
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 1122 | Line 1146 | contains  
1146      endif
1147  
1148   !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1149 <    
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 1136 | 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 = 0.0d0
1168 +    tyi = 0.0d0
1169 +    tzi = 0.0d0
1170  
1171 < !!$    write(*,*) 'sw = ', sw
1172 < !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1173 < !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1147 < !!$
1148 <    txi =  (dvduzi - dvduyi) * sw
1149 <    tyi =  (dvduxi - dvduzi) * sw
1150 <    tzi =  (dvduyi - dvduxi) * sw
1171 >    txj = 0.0d0
1172 >    tyj = 0.0d0
1173 >    tzj = 0.0d0
1174  
1175 <    txj = (dvduzj - dvduyj) * sw
1176 <    tyj = (dvduxj - dvduzj) * sw
1177 <    tzj = (dvduyj - dvduxj) * sw
1175 >    txi = (dvduyi - dvduzi) * sw
1176 >    tyi = (dvduzi - dvduxi) * sw
1177 >    tzi = (dvduxi - dvduyi) * sw
1178  
1179 < !!$    txi = -dvduxi * sw
1180 < !!$    tyi = -dvduyi * sw
1181 < !!$    tzi = -dvduzi * sw
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 <
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 1173 | 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:
1195 –    
1196 –    fxi = dvdxi * sw
1197 –    fyi = dvdyi * sw
1198 –    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 1213 | 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 1222 | 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 = 0.5_dp * (fxii + fxji)
1260 <    fyradial = 0.5_dp * (fyii + fyji)
1261 <    fzradial = 0.5_dp * (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 1243 | 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 1256 | 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  
1297    end subroutine do_shape_pair
1298 <    
1298 >
1299    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1300  
1301      ! Purpose: Compute the associated Legendre functions
# Line 1284 | Line 1313 | contains  
1313      !
1314      ! The original Fortran77 codes can be found here:
1315      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1316 <    
1316 >
1317      real (kind=dp), intent(in) :: x
1318      integer, intent(in) :: l, m, lmax
1319      real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
# Line 1301 | Line 1330 | contains  
1330  
1331      ! start with 0,0:
1332      PLM(0,0) = 1.0D0
1333 <  
1333 >
1334      ! x = +/- 1 functions are easy:
1335      IF (abs(X).EQ.1.0D0) THEN
1336         DO I = 1, m
# Line 1356 | Line 1385 | contains  
1385  
1386  
1387    subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1388 <  
1388 >
1389      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1390      !          or Ln(x) or Hn(x), and their derivatives
1391      ! Input :  function_type --- Function code
# Line 1375 | Line 1404 | contains  
1404      !
1405      ! The original Fortran77 codes can be found here:
1406      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1407 <  
1407 >
1408      real(kind=8), intent(in) :: x
1409      integer, intent(in):: m, mmax
1410      integer, intent(in):: function_type
1411      real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1412 <  
1412 >
1413      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1414      integer :: k
1415  
# Line 1426 | Line 1455 | contains  
1455  
1456  
1457      RETURN
1458 <    
1458 >
1459    end subroutine Orthogonal_Polynomial
1460 <  
1460 >
1461 >  subroutine deallocateShapes(this)
1462 >    type(Shape), pointer :: this
1463 >
1464 >    if (associated( this%ContactFuncLValue)) then
1465 >       deallocate(this%ContactFuncLValue)
1466 >       this%ContactFuncLValue => null()
1467 >    end if
1468 >
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 >    if (associated( this%ContactFuncCoefficient)) then
1479 >       deallocate(this%ContactFuncCoefficient)
1480 >       this%ContactFuncCoefficient => null()
1481 >    end if
1482 >
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 >    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