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 2188 by chuckv, Wed Apr 13 20:17:06 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 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 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 144 | 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 189 | 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 227 | 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 331 | 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 345 | Line 345 | contains  
345         status = -1
346         return
347      end if
348 <    
348 >
349      nAtypes = getSize(atypes)
350  
351      if (nAtypes == 0) then
# Line 355 | 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 371 | 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 453 | 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 489 | Line 489 | contains  
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 503 | 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 545 | 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
# Line 601 | Line 601 | contains  
601            dspidux = - (yi * xi2) / proji3
602            dspiduy = yi / proji - (yi2 * yi) / proji3
603         endif
604 <      
604 >
605         cpi = xi / proji
606         dcpidz = 0.0d0
607         dcpiduz = 0.0d0
608 <      
608 >
609         spi = yi / proji
610         dspidz = 0.0d0
611         dspiduz = 0.0d0
# Line 618 | Line 618 | contains  
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 673 | Line 673 | contains  
673                 Phunc * dlm_i(m,l) * dctidy
674            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
675                 Phunc * dlm_i(m,l) * dctidz
676 <          
676 >
677            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
678                 Phunc * dlm_i(m,l) * dctidux
679            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
# Line 688 | 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 708 | Line 708 | contains  
708            endif
709  
710            s_i = s_i + plm_i(m,l)*Phunc
711 <          
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 <          
718 >
719            dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
720                 Phunc * dlm_i(m,l) * dctidux
721            dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 724 | Line 724 | contains  
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 750 | Line 750 | contains  
750            endif
751  
752            eps_i = eps_i + plm_i(m,l)*Phunc
753 <          
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 <          
760 >
761            depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
762                 Phunc * dlm_i(m,l) * dctidux
763            depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
# Line 768 | Line 768 | contains  
768         end do
769  
770      endif
771 –      
772 –       ! 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 794 | 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 807 | 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
819         ctj = zj / rij
820 <      
820 >
821         if (ctj .gt. 1.0_dp) ctj = 1.0_dp
822         if (ctj .lt. -1.0_dp) ctj = -1.0_dp
823  
# Line 827 | Line 827 | contains  
827         dctjdux = - (zi * xj2) / r3
828         dctjduy = - (zj * yj2) / r3
829         dctjduz = zj / rij - (zj2 * zj) / r3
830 <      
830 >
831         ! this is an attempt to try to truncate the singularity when
832         ! sin(theta) is near 0.0:
833  
# Line 858 | Line 858 | contains  
858         cpj = xj / projj
859         dcpjdz = 0.0d0
860         dcpjduz = 0.0d0
861 <      
861 >
862         spj = yj / projj
863         dspjdz = 0.0d0
864         dspjduz = 0.0d0
# Line 869 | Line 869 | contains  
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 920 | 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
923 <
923 >
924            sigma_j = sigma_j + plm_j(m,l)*Phunc
925 <          
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 <          
932 >
933            dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
934                 Phunc * dlm_j(m,l) * dctjdux
935            dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
# Line 964 | Line 964 | contains  
964            endif
965  
966            s_j = s_j + plm_j(m,l)*Phunc
967 <          
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 <          
974 >
975            dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
976                 Phunc * dlm_j(m,l) * dctjdux
977            dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 1008 | Line 1008 | contains  
1008            write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1009  
1010            eps_j = eps_j + plm_j(m,l)*Phunc
1011 <          
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 <          
1018 >
1019            depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1020                 Phunc * dlm_j(m,l) * dctjdux
1021            depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
# Line 1076 | 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   !!$
# Line 1100 | 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 1122 | Line 1122 | contains  
1122      endif
1123  
1124   !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1125 <    
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 1160 | Line 1160 | contains  
1160   !!$    txj = dvduxj * sw
1161   !!$    tyj = dvduyj * sw
1162   !!$    tzj = dvduzj * sw
1163 <
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 1173 | 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 1213 | 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 1222 | Line 1222 | contains  
1222      fxij = -fxii
1223      fyij = -fyii
1224      fzij = -fzii
1225 <    
1225 >
1226      fxji = -fxjj
1227      fyji = -fyjj
1228      fzji = -fzjj
# Line 1235 | Line 1235 | contains  
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 1243 | 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 1256 | 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  
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 1284 | 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 1301 | 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 1356 | 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 1375 | 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 1426 | Line 1426 | contains  
1426  
1427  
1428      RETURN
1429 <    
1429 >
1430    end subroutine Orthogonal_Polynomial
1431 <  
1431 >
1432    subroutine deallocateShapes(this)
1433      type(Shape), pointer :: this
1434  
1435 <     if (associated( this%ContactFuncLValue)) then
1436 <        deallocate(this%ContactFuncLValue)
1437 <        this%ContactFuncLValue => null()
1438 <     end if
1439 <
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
1435 >    if (associated( this%ContactFuncLValue)) then
1436 >       deallocate(this%ContactFuncLValue)
1437 >       this%ContactFuncLValue => null()
1438 >    end if
1439  
1440 <     if (associated( this%ContactFuncCoefficient)) then
1441 <        deallocate(this%ContactFuncCoefficient)
1442 <        this%ContactFuncCoefficient => null()
1443 <     end if
1440 >    if (associated( this%ContactFuncMValue)) then
1441 >       deallocate( this%ContactFuncMValue)
1442 >       this%ContactFuncMValue => null()
1443 >    end if
1444 >    if (associated( this%ContactFunctionType)) then
1445 >       deallocate(this%ContactFunctionType)
1446 >       this%ContactFunctionType => null()
1447 >    end if
1448  
1449 <     if (associated( this%RangeFuncLValue)) then
1450 <        deallocate(this%RangeFuncLValue)
1451 <        this%RangeFuncLValue => null()
1452 <     end if
1458 <     if (associated( this%RangeFuncMValue)) then
1459 <        deallocate( this%RangeFuncMValue)
1460 <         this%RangeFuncMValue => null()
1461 <      end if
1449 >    if (associated( this%ContactFuncCoefficient)) then
1450 >       deallocate(this%ContactFuncCoefficient)
1451 >       this%ContactFuncCoefficient => null()
1452 >    end if
1453  
1454 <     if (associated( this%RangeFunctionType)) then
1455 <        deallocate( this%RangeFunctionType)
1456 <         this%RangeFunctionType => null()
1457 <      end if
1454 >    if (associated( this%RangeFuncLValue)) then
1455 >       deallocate(this%RangeFuncLValue)
1456 >       this%RangeFuncLValue => null()
1457 >    end if
1458 >    if (associated( this%RangeFuncMValue)) then
1459 >       deallocate( this%RangeFuncMValue)
1460 >       this%RangeFuncMValue => null()
1461 >    end if
1462 >
1463 >    if (associated( this%RangeFunctionType)) then
1464 >       deallocate( this%RangeFunctionType)
1465 >       this%RangeFunctionType => null()
1466 >    end if
1467      if (associated( this%RangeFuncCoefficient)) then
1468         deallocate(this%RangeFuncCoefficient)
1469         this%RangeFuncCoefficient => null()
# Line 1493 | Line 1493 | contains  
1493      integer :: i
1494      type(Shape), pointer :: thisShape
1495  
1496 < ! First walk through and kill the shape
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 <    
1501 >
1502      ! set shape map to starting values
1503      ShapeMap%n_shapes = 0
1504      ShapeMap%currentShape = 0
1505 <    
1506 <     if (associated(ShapeMap%Shapes)) then
1505 >
1506 >    if (associated(ShapeMap%Shapes)) then
1507         deallocate(ShapeMap%Shapes)
1508         ShapeMap%Shapes => null()
1509      end if
# Line 1516 | Line 1516 | contains  
1516  
1517    end subroutine destroyShapeTypes
1518  
1519 <  
1519 >
1520   end module shapes

Diff Legend

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