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 2211 by chrisfen, Thu Apr 21 14:12:19 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()
102    end type ShapeList
103 <  
103 >
104    type(ShapeList), save :: ShapeMap
105  
106    integer :: lmax
# Line 112 | Line 112 | contains  
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115 <       myATID, status)
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 166 | Line 167 | contains  
167         status = -1
168         return
169      endif
170 +    
171 +    myATID = getFirstMatchingElement(atypes, 'c_ident', c_ident)
172 + !    write(*,*) 'myATID= ', myATID, ' c_ident = ', c_ident
173  
174 <    call getElementProperty(atypes, myATID, 'c_ident', me)
175 <
172 <    ShapeMap%atidToShape(me)                         = current
173 <    ShapeMap%Shapes(current)%atid                    = me
174 >    ShapeMap%atidToShape(myATID)                     = current
175 >    ShapeMap%Shapes(current)%atid                    = myATID
176      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
177      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
178      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 189 | Line 191 | contains  
191  
192      bigL = -1
193      bigM = -1
194 <    
194 >
195      do j = 1, ShapeMap%Shapes(current)%nContactFuncs
196         if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
197            bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
# Line 227 | Line 229 | contains  
229      type(Shape), intent(inout) :: myShape
230      integer, intent(out) :: stat
231      integer :: alloc_stat
232 <
232 >
233      stat = 0
234      if (associated(myShape%contactFuncLValue)) then
235         deallocate(myShape%contactFuncLValue)
# Line 331 | Line 333 | contains  
333      return
334  
335    end subroutine allocateShape
336 <    
336 >
337    subroutine complete_Shape_FF(status)
338      integer :: status
339      integer :: i, j, l, m, lm, function_type
340      real(kind=dp) :: thisDP, sigma
341 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
341 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
342      logical :: thisProperty
343  
344      status = 0
# Line 345 | Line 347 | contains  
347         status = -1
348         return
349      end if
350 <    
350 >
351      nAtypes = getSize(atypes)
352  
353      if (nAtypes == 0) then
354         status = -1
355         return
356 <    end if
356 >    end if      
357  
358      ! atypes comes from c side
359 <    do i = 0, nAtypes
360 <      
361 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
359 >    do i = 1, nAtypes
360 >    
361 >       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
362 >       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
363 > !       write(*,*) 'is_LJ = ', thisProperty, ' for atid = ', myATID
364        
365         if (thisProperty) then
366 <          
366 >                
367            ShapeMap%currentShape = ShapeMap%currentShape + 1
368            current = ShapeMap%currentShape
369  
370 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
371 <          ShapeMap%atidToShape(thisIP) = current
368 <          ShapeMap%Shapes(current)%atid = thisIP
370 >          ShapeMap%atidToShape(myATID) = current
371 >          ShapeMap%Shapes(current)%atid = myATID
372  
373            ShapeMap%Shapes(current)%isLJ = .true.
374  
375 <          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
376 <          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
377 <          
375 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
376 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
377 >
378         endif
379 <      
379 >
380      end do
381  
382      haveShapeMap = .true.
383 <    
383 >
384 > !    do i = 1, ShapeMap%n_shapes
385 > !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
386 > !    end do
387 >
388    end subroutine complete_Shape_FF
389 <    
389 >
390    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
391         pot, A, f, t, do_pot)
392 <    
392 >
393      INTEGER, PARAMETER:: LMAX         = 64
394      INTEGER, PARAMETER:: MMAX         = 64
395  
# Line 453 | Line 460 | contains  
460      real (kind=dp) :: dsduxi, dsduyi, dsduzi
461      real (kind=dp) :: dsdxj, dsdyj, dsdzj
462      real (kind=dp) :: dsduxj, dsduyj, dsduzj
463 <    
463 >
464      real (kind=dp) :: depsdxi, depsdyi, depsdzi
465      real (kind=dp) :: depsduxi, depsduyi, depsduzi
466      real (kind=dp) :: depsdxj, depsdyj, depsdzj
# Line 489 | Line 496 | contains  
496         call handleError("calc_shape", "NO SHAPEMAP!!!!")
497         return      
498      endif
499 <    
499 >
500      !! We assume that the rotation matrices have already been calculated
501      !! and placed in the A array.
495 –
502      r3 = r2*rij
503      r5 = r3*r2
504 <    
504 >
505      drdxi = -d(1) / rij
506      drdyi = -d(2) / rij
507      drdzi = -d(3) / rij
# Line 503 | Line 509 | contains  
509      drdxj = d(1) / rij
510      drdyj = d(2) / rij
511      drdzj = d(3) / rij
512 <    
512 >
513      ! find the atom type id (atid) for each atom:
514   #ifdef IS_MPI
515      atid1 = atid_Row(atom1)
# Line 516 | Line 522 | contains  
522      ! use the atid to find the shape type (st) for each atom:
523      st1 = ShapeMap%atidToShape(atid1)
524      st2 = ShapeMap%atidToShape(atid2)
525 +    
526 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
527  
528      if (ShapeMap%Shapes(st1)%isLJ) then
529  
# Line 545 | Line 553 | contains  
553   #ifdef IS_MPI
554         ! rotate the inter-particle separation into the two different
555         ! body-fixed coordinate systems:
556 <      
556 >
557         xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
558         yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
559         zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
560 <      
560 >
561   #else
562         ! rotate the inter-particle separation into the two different
563         ! body-fixed coordinate systems:
564 <      
564 >
565         xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
566         yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
567         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
568 <      
568 >
569   #endif
570  
571         xi2 = xi*xi
# Line 601 | Line 609 | contains  
609            dspidux = - (yi * xi2) / proji3
610            dspiduy = yi / proji - (yi2 * yi) / proji3
611         endif
612 <      
612 >
613         cpi = xi / proji
614         dcpidz = 0.0d0
615         dcpiduz = 0.0d0
616 <      
616 >
617         spi = yi / proji
618         dspidz = 0.0d0
619         dspiduz = 0.0d0
# Line 618 | Line 626 | contains  
626              CHEBYSHEV_TN, tm_i, dtm_i)
627         call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
628              CHEBYSHEV_UN, um_i, dum_i)
629 <      
629 >
630         sigma_i = 0.0d0
631         s_i = 0.0d0
632         eps_i = 0.0d0
# Line 666 | Line 674 | contains  
674            endif
675  
676            sigma_i = sigma_i + plm_i(m,l)*Phunc
677 <
677 >          write(*,*) 'dsigmaidux = ', dsigmaidux
678 >          write(*,*) 'Phunc = ', Phunc
679            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
680                 Phunc * dlm_i(m,l) * dctidx
681            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
682                 Phunc * dlm_i(m,l) * dctidy
683            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
684                 Phunc * dlm_i(m,l) * dctidz
676 –          
685            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
686                 Phunc * dlm_i(m,l) * dctidux
687            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
688                 Phunc * dlm_i(m,l) * dctiduy
689            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
690                 Phunc * dlm_i(m,l) * dctiduz
691 <
691 >          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
692 >                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
693 >                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
694         end do
695  
696         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 688 | Line 698 | contains  
698            m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
699            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
700            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
701 <          
701 >
702            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
703               Phunc = coeff * tm_i(m)
704               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 708 | Line 718 | contains  
718            endif
719  
720            s_i = s_i + plm_i(m,l)*Phunc
721 <          
721 >
722            dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
723                 Phunc * dlm_i(m,l) * dctidx
724            dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
725                 Phunc * dlm_i(m,l) * dctidy
726            dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
727                 Phunc * dlm_i(m,l) * dctidz
728 <          
728 >
729            dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
730                 Phunc * dlm_i(m,l) * dctidux
731            dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 724 | Line 734 | contains  
734                 Phunc * dlm_i(m,l) * dctiduz      
735  
736         end do
737 <              
737 >
738         do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
739            l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
740            m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
741            coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
742            function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
743 <          
743 >
744            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
745               Phunc = coeff * tm_i(m)
746               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 750 | Line 760 | contains  
760            endif
761  
762            eps_i = eps_i + plm_i(m,l)*Phunc
763 <          
763 >
764            depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
765                 Phunc * dlm_i(m,l) * dctidx
766            depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
767                 Phunc * dlm_i(m,l) * dctidy
768            depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
769                 Phunc * dlm_i(m,l) * dctidz
770 <          
770 >
771            depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
772                 Phunc * dlm_i(m,l) * dctidux
773            depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 768 | Line 778 | contains  
778         end do
779  
780      endif
771 –      
772 –       ! now do j:
781  
782 +    ! now do j:
783 +
784      if (ShapeMap%Shapes(st2)%isLJ) then
785         sigma_j = ShapeMap%Shapes(st2)%sigma
786         s_j = ShapeMap%Shapes(st2)%sigma
# Line 794 | Line 804 | contains  
804         depsjduy = 0.0d0
805         depsjduz = 0.0d0
806      else
807 <      
807 >
808   #ifdef IS_MPI
809         ! rotate the inter-particle separation into the two different
810         ! body-fixed coordinate systems:
811         ! negative sign because this is the vector from j to i:
812 <      
812 >
813         xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
814         yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
815         zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
# Line 807 | Line 817 | contains  
817         ! rotate the inter-particle separation into the two different
818         ! body-fixed coordinate systems:
819         ! negative sign because this is the vector from j to i:
820 <      
820 >
821         xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
822         yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
823         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
824   #endif
825 <      
825 >
826         xj2 = xj*xj
827         yj2 = yj*yj
828         zj2 = zj*zj
829         ctj = zj / rij
830 <      
830 >
831         if (ctj .gt. 1.0_dp) ctj = 1.0_dp
832         if (ctj .lt. -1.0_dp) ctj = -1.0_dp
833  
# Line 827 | Line 837 | contains  
837         dctjdux = - (zi * xj2) / r3
838         dctjduy = - (zj * yj2) / r3
839         dctjduz = zj / rij - (zj2 * zj) / r3
840 <      
840 >
841         ! this is an attempt to try to truncate the singularity when
842         ! sin(theta) is near 0.0:
843  
# Line 858 | Line 868 | contains  
868         cpj = xj / projj
869         dcpjdz = 0.0d0
870         dcpjduz = 0.0d0
871 <      
871 >
872         spj = yj / projj
873         dspjdz = 0.0d0
874         dspjduz = 0.0d0
875  
876  
877 <       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
878 <       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
877 > !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
878 > !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
879         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
880              ShapeMap%Shapes(st2)%bigL, LMAX, &
881              plm_j, dlm_j)
882 <      
882 >
883         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
884              CHEBYSHEV_TN, tm_j, dtm_j)
885         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
886              CHEBYSHEV_UN, um_j, dum_j)
887 <      
887 >
888         sigma_j = 0.0d0
889         s_j = 0.0d0
890         eps_j = 0.0d0
# Line 920 | Line 930 | contains  
930               dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
931               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
932            endif
933 <
933 >
934            sigma_j = sigma_j + plm_j(m,l)*Phunc
935 <          
935 >
936            dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
937                 Phunc * dlm_j(m,l) * dctjdx
938            dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
939                 Phunc * dlm_j(m,l) * dctjdy
940            dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
941                 Phunc * dlm_j(m,l) * dctjdz
942 <          
942 >
943            dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
944                 Phunc * dlm_j(m,l) * dctjdux
945            dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
# Line 964 | Line 974 | contains  
974            endif
975  
976            s_j = s_j + plm_j(m,l)*Phunc
977 <          
977 >
978            dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
979                 Phunc * dlm_j(m,l) * dctjdx
980            dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
981                 Phunc * dlm_j(m,l) * dctjdy
982            dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
983                 Phunc * dlm_j(m,l) * dctjdz
984 <          
984 >
985            dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
986                 Phunc * dlm_j(m,l) * dctjdux
987            dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 1005 | Line 1015 | contains  
1015               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1016            endif
1017  
1018 <          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1018 > !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1019  
1020            eps_j = eps_j + plm_j(m,l)*Phunc
1021 <          
1021 >
1022            depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1023                 Phunc * dlm_j(m,l) * dctjdx
1024            depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1025                 Phunc * dlm_j(m,l) * dctjdy
1026            depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1027                 Phunc * dlm_j(m,l) * dctjdz
1028 <          
1028 >
1029            depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1030                 Phunc * dlm_j(m,l) * dctjdux
1031            depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 1030 | Line 1040 | contains  
1040      ! phew, now let's assemble the potential energy:
1041  
1042      sigma = 0.5*(sigma_i + sigma_j)
1043 <
1043 > !    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1044      dsigmadxi = 0.5*dsigmaidx
1045      dsigmadyi = 0.5*dsigmaidy
1046      dsigmadzi = 0.5*dsigmaidz
# Line 1062 | Line 1072 | contains  
1072      dsduzj = 0.5*dsjduz
1073  
1074      eps = sqrt(eps_i * eps_j)
1075 <
1075 >    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1076 >    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1077 > !    write(*,*) sigma_i, ' is sigma i; ', s_i, ' is s i; ', eps_i, ' is eps i'
1078      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1079      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1080      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1076 | Line 1088 | contains  
1088      depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1089      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1090      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1091 <    
1092 < !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 < !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1094 < !!$
1095 < !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1096 < !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1091 >
1092 > !    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 > !    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1094 >
1095 > !    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1096 > !    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1097   !!$
1098   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1099  
# Line 1092 | Line 1104 | contains  
1104      drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1105      drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1106      drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1107 +    write(*,*) dsduxi, ' is dsduxi; ', drduxi, ' is drduxi; ', dsigmaduxi, &
1108 +               ' is dsigmaduxi; ', dsduxi, ' is dsduxi'
1109      drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1110      drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1111      drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
# Line 1100 | Line 1114 | contains  
1114      drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1115      drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1116      drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1117 <    
1117 >
1118 > !    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1119 > !    write(*,*) 'drtdu_i = ', drtduxi, drtduyi, drtduzi
1120 >
1121      rt2 = rt*rt
1122      rt3 = rt2*rt
1123      rt5 = rt2*rt3
# Line 1122 | Line 1139 | contains  
1139      endif
1140  
1141   !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1142 <    
1142 >
1143      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1144      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1145      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1136 | Line 1153 | contains  
1153      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1154      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1155      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1156 <
1156 >    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1157      ! do the torques first since they are easy:
1158      ! remember that these are still in the body fixed axes
1159  
# Line 1145 | Line 1162 | contains  
1162   !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1163   !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1164   !!$
1165 <    txi =  (dvduzi - dvduyi) * sw
1166 <    tyi =  (dvduxi - dvduzi) * sw
1167 <    tzi =  (dvduyi - dvduxi) * sw
1165 > !    txi =  (dvduzi - dvduyi) * sw
1166 > !    tyi =  (dvduxi - dvduzi) * sw
1167 > !    tzi =  (dvduyi - dvduxi) * sw
1168  
1169 <    txj = (dvduzj - dvduyj) * sw
1170 <    tyj = (dvduxj - dvduzj) * sw
1171 <    tzj = (dvduyj - dvduxj) * sw
1169 > !    txj = (dvduzj - dvduyj) * sw
1170 > !    tyj = (dvduxj - dvduzj) * sw
1171 > !    tzj = (dvduyj - dvduxj) * sw
1172  
1173 < !!$    txi = -dvduxi * sw
1174 < !!$    tyi = -dvduyi * sw
1175 < !!$    tzi = -dvduzi * sw
1176 < !!$
1177 < !!$    txj = dvduxj * sw
1178 < !!$    tyj = dvduyj * sw
1179 < !!$    tzj = dvduzj * sw
1180 <
1173 >    txi = dvduxi * sw
1174 >    tyi = dvduyi * sw
1175 >    tzi = dvduzi * sw
1176 >
1177 >    txj = dvduxj * sw
1178 >    tyj = dvduyj * sw
1179 >    tzj = dvduzj * sw
1180 >
1181      write(*,*) 't1 = ', txi, tyi, tzi
1182      write(*,*) 't2 = ', txj, tyj, tzj
1183  
1184      ! go back to lab frame using transpose of rotation matrix:
1185 <    
1185 >
1186   #ifdef IS_MPI
1187      t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
1188           a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
# Line 1173 | Line 1190 | contains  
1190           a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
1191      t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
1192           a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
1193 <    
1193 >
1194      t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
1195           a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
1196      t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
1197 <            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1197 >         a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1198      t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
1199           a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
1200   #else
1201      t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
1202      t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
1203      t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
1204 <    
1204 >
1205      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1206      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1207      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1208   #endif
1209      ! Now, on to the forces:
1210 <    
1210 >
1211      ! first rotate the i terms back into the lab frame:
1212 <    
1212 >
1213      fxi = dvdxi * sw
1214      fyi = dvdyi * sw
1215      fzi = dvdzi * sw
# Line 1201 | Line 1218 | contains  
1218      fyj = dvdyj * sw
1219      fzj = dvdzj * sw
1220  
1221 +
1222   #ifdef IS_MPI
1223      fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1224      fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
# Line 1213 | Line 1231 | contains  
1231      fxii = a(1,atom1)*fxi + a(4,atom1)*fyi + a(7,atom1)*fzi
1232      fyii = a(2,atom1)*fxi + a(5,atom1)*fyi + a(8,atom1)*fzi
1233      fzii = a(3,atom1)*fxi + a(6,atom1)*fyi + a(9,atom1)*fzi
1234 <    
1234 >
1235      fxjj = a(1,atom2)*fxj + a(4,atom2)*fyj + a(7,atom2)*fzj
1236      fyjj = a(2,atom2)*fxj + a(5,atom2)*fyj + a(8,atom2)*fzj
1237      fzjj = a(3,atom2)*fxj + a(6,atom2)*fyj + a(9,atom2)*fzj
# Line 1222 | Line 1240 | contains  
1240      fxij = -fxii
1241      fyij = -fyii
1242      fzij = -fzii
1243 <    
1243 >
1244      fxji = -fxjj
1245      fyji = -fyjj
1246      fzji = -fzjj
# Line 1230 | Line 1248 | contains  
1248      fxradial = 0.5_dp * (fxii + fxji)
1249      fyradial = 0.5_dp * (fyii + fyji)
1250      fzradial = 0.5_dp * (fzii + fzji)
1251 <
1251 >    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1252   #ifdef IS_MPI
1253      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1254      f_Row(2,atom1) = f_Row(2,atom1) + fyradial
1255      f_Row(3,atom1) = f_Row(3,atom1) + fzradial
1256 <    
1256 >
1257      f_Col(1,atom2) = f_Col(1,atom2) - fxradial
1258      f_Col(2,atom2) = f_Col(2,atom2) - fyradial
1259      f_Col(3,atom2) = f_Col(3,atom2) - fzradial
# Line 1243 | Line 1261 | contains  
1261      f(1,atom1) = f(1,atom1) + fxradial
1262      f(2,atom1) = f(2,atom1) + fyradial
1263      f(3,atom1) = f(3,atom1) + fzradial
1264 <    
1264 >
1265      f(1,atom2) = f(1,atom2) - fxradial
1266      f(2,atom2) = f(2,atom2) - fyradial
1267      f(3,atom2) = f(3,atom2) - fzradial
# Line 1256 | Line 1274 | contains  
1274      id1 = atom1
1275      id2 = atom2
1276   #endif
1277 <    
1277 >
1278      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1279 <      
1279 >
1280         fpair(1) = fpair(1) + fxradial
1281         fpair(2) = fpair(2) + fyradial
1282         fpair(3) = fpair(3) + fzradial
1283 <      
1283 >
1284      endif
1285  
1286    end subroutine do_shape_pair
1287 <    
1287 >
1288    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1289  
1290      ! Purpose: Compute the associated Legendre functions
# Line 1284 | Line 1302 | contains  
1302      !
1303      ! The original Fortran77 codes can be found here:
1304      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1305 <    
1305 >
1306      real (kind=dp), intent(in) :: x
1307      integer, intent(in) :: l, m, lmax
1308      real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
# Line 1301 | Line 1319 | contains  
1319  
1320      ! start with 0,0:
1321      PLM(0,0) = 1.0D0
1322 <  
1322 >
1323      ! x = +/- 1 functions are easy:
1324      IF (abs(X).EQ.1.0D0) THEN
1325         DO I = 1, m
# Line 1356 | Line 1374 | contains  
1374  
1375  
1376    subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1377 <  
1377 >
1378      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1379      !          or Ln(x) or Hn(x), and their derivatives
1380      ! Input :  function_type --- Function code
# Line 1375 | Line 1393 | contains  
1393      !
1394      ! The original Fortran77 codes can be found here:
1395      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1396 <  
1396 >
1397      real(kind=8), intent(in) :: x
1398      integer, intent(in):: m, mmax
1399      integer, intent(in):: function_type
1400      real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1401 <  
1401 >
1402      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1403      integer :: k
1404  
# Line 1426 | Line 1444 | contains  
1444  
1445  
1446      RETURN
1447 <    
1447 >
1448    end subroutine Orthogonal_Polynomial
1449 <  
1449 >
1450 >  subroutine deallocateShapes(this)
1451 >    type(Shape), pointer :: this
1452 >
1453 >    if (associated( this%ContactFuncLValue)) then
1454 >       deallocate(this%ContactFuncLValue)
1455 >       this%ContactFuncLValue => null()
1456 >    end if
1457 >
1458 >    if (associated( this%ContactFuncMValue)) then
1459 >       deallocate( this%ContactFuncMValue)
1460 >       this%ContactFuncMValue => null()
1461 >    end if
1462 >    if (associated( this%ContactFunctionType)) then
1463 >       deallocate(this%ContactFunctionType)
1464 >       this%ContactFunctionType => null()
1465 >    end if
1466 >
1467 >    if (associated( this%ContactFuncCoefficient)) then
1468 >       deallocate(this%ContactFuncCoefficient)
1469 >       this%ContactFuncCoefficient => null()
1470 >    end if
1471 >
1472 >    if (associated( this%RangeFuncLValue)) then
1473 >       deallocate(this%RangeFuncLValue)
1474 >       this%RangeFuncLValue => null()
1475 >    end if
1476 >    if (associated( this%RangeFuncMValue)) then
1477 >       deallocate( this%RangeFuncMValue)
1478 >       this%RangeFuncMValue => null()
1479 >    end if
1480 >
1481 >    if (associated( this%RangeFunctionType)) then
1482 >       deallocate( this%RangeFunctionType)
1483 >       this%RangeFunctionType => null()
1484 >    end if
1485 >    if (associated( this%RangeFuncCoefficient)) then
1486 >       deallocate(this%RangeFuncCoefficient)
1487 >       this%RangeFuncCoefficient => null()
1488 >    end if
1489 >
1490 >    if (associated( this%StrengthFuncLValue)) then
1491 >       deallocate(this%StrengthFuncLValue)
1492 >       this%StrengthFuncLValue => null()
1493 >    end if
1494 >
1495 >    if (associated( this%StrengthFuncMValue )) then
1496 >       deallocate(this%StrengthFuncMValue)
1497 >       this%StrengthFuncMValue => null()
1498 >    end if
1499 >
1500 >    if(associated( this%StrengthFunctionType)) then
1501 >       deallocate(this%StrengthFunctionType)
1502 >       this%StrengthFunctionType => null()
1503 >    end if
1504 >    if (associated( this%StrengthFuncCoefficient )) then
1505 >       deallocate(this%StrengthFuncCoefficient)
1506 >       this%StrengthFuncCoefficient => null()
1507 >    end if
1508 >  end subroutine deallocateShapes
1509 >
1510 >  subroutine destroyShapeTypes
1511 >    integer :: i
1512 >    type(Shape), pointer :: thisShape
1513 >
1514 >    ! First walk through and kill the shape
1515 >    do i = 1,ShapeMap%n_shapes
1516 >       thisShape => ShapeMap%Shapes(i)
1517 >       call deallocateShapes(thisShape)
1518 >    end do
1519 >
1520 >    ! set shape map to starting values
1521 >    ShapeMap%n_shapes = 0
1522 >    ShapeMap%currentShape = 0
1523 >
1524 >    if (associated(ShapeMap%Shapes)) then
1525 >       deallocate(ShapeMap%Shapes)
1526 >       ShapeMap%Shapes => null()
1527 >    end if
1528 >
1529 >    if (associated(ShapeMap%atidToShape)) then
1530 >       deallocate(ShapeMap%atidToShape)
1531 >       ShapeMap%atidToShape => null()
1532 >    end if
1533 >
1534 >
1535 >  end subroutine destroyShapeTypes
1536 >
1537 >
1538   end module shapes

Diff Legend

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