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

Comparing trunk/OOPSE-2.0/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1688 by chrisfen, Tue Oct 26 18:03:12 2004 UTC vs.
Revision 1689 by chrisfen, Fri Oct 29 22:28:45 2004 UTC

# Line 73 | Line 73 | contains  
73         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
74         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
75         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
76 <       myAtid, status)
76 >       myATID, status)
77  
78      integer :: nContactFuncs
79      integer :: nRangeFuncs
80      integer :: nStrengthFuncs
81      integer :: shape_ident
82      integer :: status
83 <    integer :: myAtid
83 >    integer :: myATID
84      integer :: bigL
85      integer :: bigM
86      integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
# Line 115 | Line 115 | contains  
115        
116         ntypes = getSize(atypes)
117        
118 <       allocate(ShapeMap%atidToShape(ntypes))
118 >       allocate(ShapeMap%atidToShape(0:ntypes))
119      end if
120      
121      ShapeMap%currentShape = ShapeMap%currentShape + 1
# Line 128 | Line 128 | contains  
128         return
129      endif
130  
131 <    call getElementProperty(atypes, myAtid, "c_ident", me)
131 >    call getElementProperty(atypes, myATID, 'c_ident', me)
132 >
133      ShapeMap%atidToShape(me)                         = current
134      ShapeMap%Shapes(current)%atid                    = me
135      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
# Line 188 | Line 189 | contains  
189      integer, intent(out) :: stat
190      integer :: alloc_stat
191  
192 +    stat = 0
193      if (associated(myShape%contactFuncLValue)) then
194         deallocate(myShape%contactFuncLValue)
195      endif
# Line 253 | Line 255 | contains  
255         stat = -1
256         return
257      endif
258 <    
258 >
259      if (associated(myShape%strengthFuncLValue)) then
260         deallocate(myShape%strengthFuncLValue)
261      endif
# Line 286 | Line 288 | contains  
288         stat = -1
289         return
290      endif
291 +
292 +    return
293  
294    end subroutine allocateShape
295      
# Line 310 | Line 314 | contains  
314         return
315      end if
316  
317 <    do i = 1, nAtypes
317 >    ! atypes comes from c side
318 >    do i = 0, nAtypes
319        
320         call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
321        
# Line 339 | Line 344 | contains  
344    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
345         pot, A, f, t, do_pot)
346      
347 +    INTEGER, PARAMETER:: LMAX         = 64
348 +    INTEGER, PARAMETER:: MMAX         = 64
349 +
350      integer, intent(in) :: atom1, atom2
351      real (kind=dp), intent(inout) :: rij, r2
352      real (kind=dp), dimension(3), intent(in) :: d
# Line 430 | Line 438 | contains  
438      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
439      real (kind=dp) :: fxradial, fyradial, fzradial
440  
441 +    real (kind=dp) :: plm_i(LMAX,MMAX), dlm_i(LMAX,MMAX)
442 +    real (kind=dp) :: plm_j(LMAX,MMAX), dlm_j(LMAX,MMAX)
443 +    real (kind=dp) :: tm_i(MMAX), dtm_i(MMAX), um_i(MMAX), dum_i(MMAX)
444 +    real (kind=dp) :: tm_j(MMAX), dtm_j(MMAX), um_j(MMAX), dum_j(MMAX)
445 +
446      if (.not.haveShapeMap) then
447         call handleError("calc_shape", "NO SHAPEMAP!!!!")
448         return      
# Line 437 | Line 450 | contains  
450      
451      !! We assume that the rotation matrices have already been calculated
452      !! and placed in the A array.
453 <        
453 >
454      r3 = r2*rij
455      r5 = r3*r2
456      
# Line 459 | Line 472 | contains  
472   #endif
473  
474      ! use the atid to find the shape type (st) for each atom:
462 –
475      st1 = ShapeMap%atidToShape(atid1)
476      st2 = ShapeMap%atidToShape(atid2)
477 <    
477 >
478      if (ShapeMap%Shapes(st1)%isLJ) then
479 +
480         sigma_i = ShapeMap%Shapes(st1)%sigma
481         s_i = ShapeMap%Shapes(st1)%sigma
482         eps_i = ShapeMap%Shapes(st1)%epsilon
# Line 504 | Line 517 | contains  
517         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
518        
519   #endif
520 <      
520 >
521         xi2 = xi*xi
522         yi2 = yi*yi
523         zi2 = zi*zi
# Line 513 | Line 526 | contains  
526         proji3 = proji*proji*proji
527        
528         cti = zi / rij
529 +
530         dctidx = - zi * xi / r3
531         dctidy = - zi * yi / r3
532         dctidz = 1.0d0 / rij - zi2 / r3
# Line 537 | Line 551 | contains  
551         dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
552  
553         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
554 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
554 >            ShapeMap%Shapes(st1)%bigM, LMAX, &
555 >            plm_i, dlm_i)
556  
557 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
557 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
558              CHEBYSHEV_TN, tm_i, dtm_i)
559 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
559 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
560              CHEBYSHEV_UN, um_i, dum_i)
561        
562         sigma_i = 0.0d0
# Line 573 | Line 588 | contains  
588            function_type = ShapeMap%Shapes(st1)%ContactFunctionType(lm)
589  
590            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
591 +           !  write(*,*) tm_i(m), ' is tm_i'
592               Phunc = coeff * tm_i(m)
593               dPhuncdX = coeff * dtm_i(m) * dcpidx
594               dPhuncdY = coeff * dtm_i(m) * dcpidy
# Line 591 | Line 607 | contains  
607            endif
608  
609            sigma_i = sigma_i + plm_i(l,m)*Phunc
610 <          
610 >          write(*,*) plm_i(l,m), l, m
611            dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
612                 Phunc * dlm_i(l,m) * dctidx
613            dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
# Line 673 | Line 689 | contains  
689               dPhuncdUy = coeff*(spi * dum_i(m-1)*dcpiduy + dspiduy *um_i(m-1))
690               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
691            endif
692 <
692 >          !write(*,*) eps_i, plm_i(l,m), Phunc
693            eps_i = eps_i + plm_i(l,m)*Phunc
694            
695            depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
# Line 770 | Line 786 | contains  
786         dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
787        
788         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
789 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
790 <      
791 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
789 >            ShapeMap%Shapes(st2)%bigM, LMAX, &
790 >            plm_j, dlm_j)
791 >      
792 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
793              CHEBYSHEV_TN, tm_j, dtm_j)
794 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
794 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
795              CHEBYSHEV_UN, um_j, dum_j)
796        
797         sigma_j = 0.0d0
# Line 960 | Line 977 | contains  
977      dsduxj = 0.5*dsjdux
978      dsduyj = 0.5*dsjduy
979      dsduzj = 0.5*dsjduz
980 <
980 >    !write(*,*) eps_i, eps_j
981      eps = sqrt(eps_i * eps_j)
982  
983      depsdxi = eps_j * depsidx / (2.0d0 * eps)
# Line 1138 | Line 1155 | contains  
1155      
1156    end subroutine do_shape_pair
1157      
1158 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1159 <    
1158 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1159 >
1160      ! Purpose: Compute the associated Legendre functions
1161      !          Plm(x) and their derivatives Plm'(x)
1162      ! Input :  x  --- Argument of Plm(x)
# Line 1156 | Line 1173 | contains  
1173      ! The original Fortran77 codes can be found here:
1174      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1175      
1176 <    real (kind=8), intent(in) :: x
1176 >    real (kind=dp), intent(in) :: x
1177      integer, intent(in) :: l, m, lmax
1178 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1178 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1179      integer :: i, j, ls
1180 <    real (kind=8) :: xq, xs
1180 >    real (kind=dp) :: xq, xs
1181  
1182      ! zero out both arrays:
1183      DO I = 0, m
1184         DO J = 0, l
1185 <          PLM(J,I) = 0.0D0
1186 <          DLM(J,I) = 0.0D0
1185 >          PLM(J,I) = 0.0_dp
1186 >          DLM(J,I) = 0.0_dp
1187         end DO
1188      end DO
1189  
# Line 1199 | Line 1216 | contains  
1216      DO I = 1, l
1217         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1218      enddo
1219 <    
1219 >
1220      DO I = 0, l
1221         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1222      enddo
1223 <    
1223 >
1224      DO I = 0, l
1225         DO J = I+2, m
1226            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1227                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1228         end DO
1229      end DO
1230 <  
1230 >
1231      DLM(0, 0)=0.0D0
1215 –    
1232      DO J = 1, m
1233         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1234      end DO
1235 <    
1235 >
1236      DO I = 1, l
1237         DO J = I, m
1238            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1239         end DO
1240      end DO
1241 <    
1241 >
1242      RETURN
1243    END SUBROUTINE Associated_Legendre
1244  
1245  
1246 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1246 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1247    
1248      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1249      !          or Ln(x) or Hn(x), and their derivatives
# Line 1249 | Line 1265 | contains  
1265      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1266    
1267      real(kind=8), intent(in) :: x
1268 <    integer, intent(in):: m
1268 >    integer, intent(in):: m, mmax
1269      integer, intent(in):: function_type
1270 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1270 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1271    
1272      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1273      integer :: k
# Line 1306 | Line 1322 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1322       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1323       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1324       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1325 <     myAtid, status)
1325 >     myATID, status)
1326  
1327    use definitions
1328    use shapes, only: newShapeType
# Line 1315 | Line 1331 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1331    integer :: nRangeFuncs
1332    integer :: nStrengthFuncs
1333    integer :: status
1334 <  integer :: myAtid
1334 >  integer :: myATID
1335    
1336    integer, dimension(nContactFuncs) :: ContactFuncLValue          
1337    integer, dimension(nContactFuncs) :: ContactFuncMValue          
# Line 1335 | Line 1351 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1351         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1352         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1353         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1354 <       myAtid, status)
1354 >       myATID, status)
1355  
1356    return
1357   end subroutine makeShape

Diff Legend

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