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

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1608 by gezelter, Wed Oct 20 04:02:48 2004 UTC vs.
Revision 1698 by chrisfen, Tue Nov 2 15:28:25 2004 UTC

# Line 6 | Line 6 | module shapes
6    use vector_class
7    use simulation
8    use status
9 +  use lj
10   #ifdef IS_MPI
11    use mpiSimulation
12   #endif
# Line 24 | Line 25 | module shapes
25  
26    public :: do_shape_pair
27    public :: newShapeType
28 +  public :: complete_Shape_FF
29  
30  
31    type, private :: Shape
# Line 60 | Line 62 | module shapes
62    type(ShapeList), save :: ShapeMap
63  
64    integer :: lmax
63 –  real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
64 –  real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
65 –  real (kind=dp), allocatable, dimension(:) :: tm_j, dtm_j, um_j, dum_j
65  
66   contains  
67  
# Line 71 | Line 70 | contains  
70         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
71         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
72         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
73 <       myAtid, status)
73 >       myATID, status)
74  
75      integer :: nContactFuncs
76      integer :: nRangeFuncs
77      integer :: nStrengthFuncs
78      integer :: shape_ident
79      integer :: status
80 <    integer :: myAtid
80 >    integer :: myATID
81      integer :: bigL
82      integer :: bigM
83      integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
# Line 104 | Line 103 | contains  
103         call getMatchingElementList(atypes, "is_Shape", .true., &
104              nShapeTypes, MatchList)
105        
106 <       call getMatchingElementList(atypes, "is_LJ", .true., &
106 >       call getMatchingElementList(atypes, "is_LennardJones", .true., &
107              nLJTypes, MatchList)
108        
109         ShapeMap%n_shapes = nShapeTypes + nLJTypes
# Line 113 | Line 112 | contains  
112        
113         ntypes = getSize(atypes)
114        
115 <       allocate(ShapeMap%atidToShape(ntypes))
115 >       allocate(ShapeMap%atidToShape(0:ntypes))
116      end if
117      
118      ShapeMap%currentShape = ShapeMap%currentShape + 1
# Line 126 | Line 125 | contains  
125         return
126      endif
127  
128 <    call getElementProperty(atypes, myAtid, "c_ident", me)
128 >    call getElementProperty(atypes, myATID, 'c_ident', me)
129 >
130      ShapeMap%atidToShape(me)                         = current
131      ShapeMap%Shapes(current)%atid                    = me
132      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
# Line 186 | Line 186 | contains  
186      integer, intent(out) :: stat
187      integer :: alloc_stat
188  
189 +    stat = 0
190      if (associated(myShape%contactFuncLValue)) then
191         deallocate(myShape%contactFuncLValue)
192      endif
# Line 251 | Line 252 | contains  
252         stat = -1
253         return
254      endif
255 <    
255 >
256      if (associated(myShape%strengthFuncLValue)) then
257         deallocate(myShape%strengthFuncLValue)
258      endif
# Line 285 | Line 286 | contains  
286         return
287      endif
288  
289 +    return
290 +
291    end subroutine allocateShape
292      
293 <  subroutine init_Shape_FF(status)
293 >  subroutine complete_Shape_FF(status)
294      integer :: status
295      integer :: i, j, l, m, lm, function_type
296 <    real(kind=dp) :: bigSigma, myBigSigma, thisSigma, coeff, Phunc, spi
294 <    real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP
296 >    real(kind=dp) :: thisDP, sigma
297      integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
298      logical :: thisProperty
299  
298 –    Pi = 4.0d0 * datan(1.0d0)
299 –
300      status = 0
301      if (ShapeMap%currentShape == 0) then
302         call handleError("init_Shape_FF", "No members in ShapeMap")
303         status = -1
304         return
305      end if
306 <
307 <    bigSigma = 0.0d0
308 <    do i = 1, ShapeMap%currentShape
309 <
310 <       ! Scan over theta and phi to
311 <       ! find the largest contact in any direction....
312 <
313 <       myBigSigma = 0.0d0
314 <
315 <       do iTheta = 0, nSteps
316 <          theta = (Pi/2.0d0)*(dble(iTheta)/dble(nSteps))
317 <          costheta = cos(theta)
318 <
319 <          call Associated_Legendre(costheta, ShapeMap%Shapes(i)%bigL, &
320 <               ShapeMap%Shapes(i)%bigM, lmax, plm_i, dlm_i)
321 <                  
322 <          do iPhi = 0, nSteps
323 <             phi = -Pi  + 2.0d0 * Pi * (dble(iPhi)/dble(nSteps))
324 <             cpi = cos(phi)
325 <             spi = sin(phi)
326 <            
327 <             call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
328 <                  CHEBYSHEV_TN, tm_i, dtm_i)
329 <             call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
330 <                  CHEBYSHEV_UN, um_i, dum_i)
331 <
332 <             thisSigma = 0.0d0
333 <
334 <             do lm = 1, ShapeMap%Shapes(i)%nContactFuncs                
335 <
336 <                l = ShapeMap%Shapes(i)%ContactFuncLValue(lm)
337 <                m = ShapeMap%Shapes(i)%ContactFuncMValue(lm)
338 <                coeff = ShapeMap%Shapes(i)%ContactFuncCoefficient(lm)
339 <                function_type = ShapeMap%Shapes(i)%ContactFunctionType(lm)
340 <
341 <                if ((function_type .eq. SH_COS).or.(m.eq.0)) then
342 <                   Phunc = coeff * tm_i(m)
343 <                else
344 <                   Phunc = coeff * spi * um_i(m-1)
345 <                endif
346 <                
347 <                thisSigma = thisSigma + plm_i(l,m)*Phunc
348 <             enddo
349 <
350 <             if (thisSigma.gt.myBigSigma) myBigSigma = thisSigma
351 <          enddo
352 <       enddo
353 <
354 <       if (myBigSigma.gt.bigSigma) bigSigma = myBigSigma
355 <    enddo
356 <
306 >    
307      nAtypes = getSize(atypes)
308  
309      if (nAtypes == 0) then
# Line 361 | Line 311 | contains  
311         return
312      end if
313  
314 <    do i = 1, nAtypes
314 >    ! atypes comes from c side
315 >    do i = 0, nAtypes
316        
317 <       call getElementProperty(atypes, i, "is_LJ", thisProperty)
317 >       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
318        
319         if (thisProperty) then
320            
# Line 376 | Line 327 | contains  
327  
328            ShapeMap%Shapes(current)%isLJ = .true.
329  
330 <          call getElementProperty(atypes, i, "lj_epsilon", thisDP)
331 <          ShapeMap%Shapes(current)%epsilon = thisDP
381 <
382 <          call getElementProperty(atypes, i, "lj_sigma",   thisDP)
383 <          ShapeMap%Shapes(current)%sigma = thisDP          
384 <          if (thisDP .gt. bigSigma) bigSigma = thisDP
330 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
331 >          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
332            
333         endif
334        
# Line 389 | Line 336 | contains  
336  
337      haveShapeMap = .true.
338      
339 <  end subroutine init_Shape_FF
339 >  end subroutine complete_Shape_FF
340      
341    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
342         pot, A, f, t, do_pot)
343      
344 +    INTEGER, PARAMETER:: LMAX         = 64
345 +    INTEGER, PARAMETER:: MMAX         = 64
346 +
347      integer, intent(in) :: atom1, atom2
348      real (kind=dp), intent(inout) :: rij, r2
349      real (kind=dp), dimension(3), intent(in) :: d
# Line 485 | Line 435 | contains  
435      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
436      real (kind=dp) :: fxradial, fyradial, fzradial
437  
438 +    real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
439 +    real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
440 +    real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
441 +    real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
442 +
443      if (.not.haveShapeMap) then
444         call handleError("calc_shape", "NO SHAPEMAP!!!!")
445         return      
# Line 492 | Line 447 | contains  
447      
448      !! We assume that the rotation matrices have already been calculated
449      !! and placed in the A array.
450 <        
450 >
451      r3 = r2*rij
452      r5 = r3*r2
453      
# Line 514 | Line 469 | contains  
469   #endif
470  
471      ! use the atid to find the shape type (st) for each atom:
517 –
472      st1 = ShapeMap%atidToShape(atid1)
473      st2 = ShapeMap%atidToShape(atid2)
474 <    
474 >
475      if (ShapeMap%Shapes(st1)%isLJ) then
476 +
477         sigma_i = ShapeMap%Shapes(st1)%sigma
478         s_i = ShapeMap%Shapes(st1)%sigma
479         eps_i = ShapeMap%Shapes(st1)%epsilon
# Line 559 | Line 514 | contains  
514         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
515        
516   #endif
517 <      
517 >
518         xi2 = xi*xi
519         yi2 = yi*yi
520         zi2 = zi*zi
# Line 568 | Line 523 | contains  
523         proji3 = proji*proji*proji
524        
525         cti = zi / rij
526 +
527         dctidx = - zi * xi / r3
528         dctidy = - zi * yi / r3
529         dctidz = 1.0d0 / rij - zi2 / r3
# Line 591 | Line 547 | contains  
547         dspiduy = xi * yi * zi / proji3
548         dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
549  
550 <       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
551 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
550 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
551 >            ShapeMap%Shapes(st1)%bigL, LMAX, &
552 >            plm_i, dlm_i)
553  
554 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
554 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
555              CHEBYSHEV_TN, tm_i, dtm_i)
556 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
556 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
557              CHEBYSHEV_UN, um_i, dum_i)
558        
559         sigma_i = 0.0d0
# Line 645 | Line 602 | contains  
602               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
603            endif
604  
605 <          sigma_i = sigma_i + plm_i(l,m)*Phunc
605 >          sigma_i = sigma_i + plm_i(m,l)*Phunc
606 >
607 >          dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
608 >               Phunc * dlm_i(m,l) * dctidx
609 >          dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
610 >               Phunc * dlm_i(m,l) * dctidy
611 >          dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
612 >               Phunc * dlm_i(m,l) * dctidz
613            
614 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
615 <               Phunc * dlm_i(l,m) * dctidx
616 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
617 <               Phunc * dlm_i(l,m) * dctidy
618 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
619 <               Phunc * dlm_i(l,m) * dctidz
656 <          
657 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
658 <               Phunc * dlm_i(l,m) * dctidux
659 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
660 <               Phunc * dlm_i(l,m) * dctiduy
661 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
662 <               Phunc * dlm_i(l,m) * dctiduz
614 >          dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
615 >               Phunc * dlm_i(m,l) * dctidux
616 >          dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
617 >               Phunc * dlm_i(m,l) * dctiduy
618 >          dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
619 >               Phunc * dlm_i(m,l) * dctiduz
620  
621         end do
622  
# Line 687 | Line 644 | contains  
644               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
645            endif
646  
647 <          s_i = s_i + plm_i(l,m)*Phunc
647 >          s_i = s_i + plm_i(m,l)*Phunc
648            
649 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
650 <               Phunc * dlm_i(l,m) * dctidx
651 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
652 <               Phunc * dlm_i(l,m) * dctidy
653 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
654 <               Phunc * dlm_i(l,m) * dctidz
649 >          dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
650 >               Phunc * dlm_i(m,l) * dctidx
651 >          dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
652 >               Phunc * dlm_i(m,l) * dctidy
653 >          dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
654 >               Phunc * dlm_i(m,l) * dctidz
655            
656 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
657 <               Phunc * dlm_i(l,m) * dctidux
658 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
659 <               Phunc * dlm_i(l,m) * dctiduy
660 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
661 <               Phunc * dlm_i(l,m) * dctiduz      
656 >          dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
657 >               Phunc * dlm_i(m,l) * dctidux
658 >          dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
659 >               Phunc * dlm_i(m,l) * dctiduy
660 >          dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
661 >               Phunc * dlm_i(m,l) * dctiduz      
662  
663         end do
664                
# Line 729 | Line 686 | contains  
686               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
687            endif
688  
689 <          eps_i = eps_i + plm_i(l,m)*Phunc
689 >          eps_i = eps_i + plm_i(m,l)*Phunc
690            
691 <          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
692 <               Phunc * dlm_i(l,m) * dctidx
693 <          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
694 <               Phunc * dlm_i(l,m) * dctidy
695 <          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
696 <               Phunc * dlm_i(l,m) * dctidz
691 >          depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
692 >               Phunc * dlm_i(m,l) * dctidx
693 >          depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
694 >               Phunc * dlm_i(m,l) * dctidy
695 >          depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
696 >               Phunc * dlm_i(m,l) * dctidz
697            
698 <          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
699 <               Phunc * dlm_i(l,m) * dctidux
700 <          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
701 <               Phunc * dlm_i(l,m) * dctiduy
702 <          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
703 <               Phunc * dlm_i(l,m) * dctiduz      
698 >          depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
699 >               Phunc * dlm_i(m,l) * dctidux
700 >          depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
701 >               Phunc * dlm_i(m,l) * dctiduy
702 >          depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
703 >               Phunc * dlm_i(m,l) * dctiduz      
704  
705         end do
706  
# Line 824 | Line 781 | contains  
781         dspjduy = xj * yj * zj / projj3
782         dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
783        
784 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
785 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
784 >       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
785 >            ShapeMap%Shapes(st2)%bigL, LMAX, &
786 >            plm_j, dlm_j)
787        
788 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
788 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
789              CHEBYSHEV_TN, tm_j, dtm_j)
790 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
790 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
791              CHEBYSHEV_UN, um_j, dum_j)
792        
793         sigma_j = 0.0d0
# Line 878 | Line 836 | contains  
836               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
837            endif
838  
839 <          sigma_j = sigma_j + plm_j(l,m)*Phunc
839 >          sigma_j = sigma_j + plm_j(m,l)*Phunc
840            
841 <          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
842 <               Phunc * dlm_j(l,m) * dctjdx
843 <          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
844 <               Phunc * dlm_j(l,m) * dctjdy
845 <          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
846 <               Phunc * dlm_j(l,m) * dctjdz
841 >          dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
842 >               Phunc * dlm_j(m,l) * dctjdx
843 >          dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
844 >               Phunc * dlm_j(m,l) * dctjdy
845 >          dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
846 >               Phunc * dlm_j(m,l) * dctjdz
847            
848 <          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
849 <               Phunc * dlm_j(l,m) * dctjdux
850 <          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
851 <               Phunc * dlm_j(l,m) * dctjduy
852 <          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
853 <               Phunc * dlm_j(l,m) * dctjduz
848 >          dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
849 >               Phunc * dlm_j(m,l) * dctjdux
850 >          dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
851 >               Phunc * dlm_j(m,l) * dctjduy
852 >          dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
853 >               Phunc * dlm_j(m,l) * dctjduz
854  
855         end do
856  
# Line 920 | Line 878 | contains  
878               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
879            endif
880  
881 <          s_j = s_j + plm_j(l,m)*Phunc
881 >          s_j = s_j + plm_j(m,l)*Phunc
882            
883 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
884 <               Phunc * dlm_j(l,m) * dctjdx
885 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
886 <               Phunc * dlm_j(l,m) * dctjdy
887 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
888 <               Phunc * dlm_j(l,m) * dctjdz
883 >          dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
884 >               Phunc * dlm_j(m,l) * dctjdx
885 >          dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
886 >               Phunc * dlm_j(m,l) * dctjdy
887 >          dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
888 >               Phunc * dlm_j(m,l) * dctjdz
889            
890 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
891 <               Phunc * dlm_j(l,m) * dctjdux
892 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
893 <               Phunc * dlm_j(l,m) * dctjduy
894 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
895 <               Phunc * dlm_j(l,m) * dctjduz
890 >          dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
891 >               Phunc * dlm_j(m,l) * dctjdux
892 >          dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
893 >               Phunc * dlm_j(m,l) * dctjduy
894 >          dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
895 >               Phunc * dlm_j(m,l) * dctjduz
896  
897         end do
898  
# Line 962 | Line 920 | contains  
920               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
921            endif
922  
923 <          eps_j = eps_j + plm_j(l,m)*Phunc
923 >          eps_j = eps_j + plm_j(m,l)*Phunc
924            
925 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
926 <               Phunc * dlm_j(l,m) * dctjdx
927 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
928 <               Phunc * dlm_j(l,m) * dctjdy
929 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
930 <               Phunc * dlm_j(l,m) * dctjdz
925 >          depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
926 >               Phunc * dlm_j(m,l) * dctjdx
927 >          depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
928 >               Phunc * dlm_j(m,l) * dctjdy
929 >          depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
930 >               Phunc * dlm_j(m,l) * dctjdz
931            
932 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
933 <               Phunc * dlm_j(l,m) * dctjdux
934 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
935 <               Phunc * dlm_j(l,m) * dctjduy
936 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
937 <               Phunc * dlm_j(l,m) * dctjduz
932 >          depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
933 >               Phunc * dlm_j(m,l) * dctjdux
934 >          depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
935 >               Phunc * dlm_j(m,l) * dctjduy
936 >          depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
937 >               Phunc * dlm_j(m,l) * dctjduz
938  
939         end do
940  
# Line 1193 | Line 1151 | contains  
1151      
1152    end subroutine do_shape_pair
1153      
1154 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1155 <    
1154 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1155 >
1156      ! Purpose: Compute the associated Legendre functions
1157      !          Plm(x) and their derivatives Plm'(x)
1158      ! Input :  x  --- Argument of Plm(x)
# Line 1211 | Line 1169 | contains  
1169      ! The original Fortran77 codes can be found here:
1170      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1171      
1172 <    real (kind=8), intent(in) :: x
1172 >    real (kind=dp), intent(in) :: x
1173      integer, intent(in) :: l, m, lmax
1174 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1174 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1175      integer :: i, j, ls
1176 <    real (kind=8) :: xq, xs
1176 >    real (kind=dp) :: xq, xs
1177  
1178      ! zero out both arrays:
1179      DO I = 0, m
1180         DO J = 0, l
1181 <          PLM(J,I) = 0.0D0
1182 <          DLM(J,I) = 0.0D0
1181 >          PLM(J,I) = 0.0_dp
1182 >          DLM(J,I) = 0.0_dp
1183         end DO
1184      end DO
1185  
# Line 1254 | Line 1212 | contains  
1212      DO I = 1, l
1213         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1214      enddo
1215 <    
1215 >
1216      DO I = 0, l
1217         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1218      enddo
1219 <    
1219 >
1220      DO I = 0, l
1221         DO J = I+2, m
1222            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1223                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1224         end DO
1225      end DO
1226 <  
1226 >
1227      DLM(0, 0)=0.0D0
1270 –    
1228      DO J = 1, m
1229         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1230      end DO
1231 <    
1231 >
1232      DO I = 1, l
1233         DO J = I, m
1234            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1235         end DO
1236      end DO
1237 <    
1237 >
1238      RETURN
1239    END SUBROUTINE Associated_Legendre
1240  
1241  
1242 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1242 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1243    
1244      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1245      !          or Ln(x) or Hn(x), and their derivatives
# Line 1304 | Line 1261 | contains  
1261      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1262    
1263      real(kind=8), intent(in) :: x
1264 <    integer, intent(in):: m
1264 >    integer, intent(in):: m, mmax
1265      integer, intent(in):: function_type
1266 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1266 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1267    
1268      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1269      integer :: k
# Line 1350 | Line 1307 | contains  
1307         DY0 = DY1
1308         DY1 = DYN
1309      end DO
1310 +
1311 +
1312      RETURN
1313      
1314    end subroutine Orthogonal_Polynomial
# Line 1361 | Line 1320 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1320       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1321       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1322       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1323 <     myAtid, status)
1323 >     myATID, status)
1324  
1325    use definitions
1326    use shapes, only: newShapeType
# Line 1370 | Line 1329 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1329    integer :: nRangeFuncs
1330    integer :: nStrengthFuncs
1331    integer :: status
1332 <  integer :: myAtid
1332 >  integer :: myATID
1333    
1334    integer, dimension(nContactFuncs) :: ContactFuncLValue          
1335    integer, dimension(nContactFuncs) :: ContactFuncMValue          
# Line 1390 | Line 1349 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1349         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1350         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1351         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1352 <       myAtid, status)
1352 >       myATID, status)
1353  
1354    return
1355   end subroutine makeShape
1356 +
1357 + subroutine completeShapeFF(status)
1358 +
1359 +  use shapes, only: complete_Shape_FF
1360 +
1361 +  integer, intent(out)  :: status
1362 +  integer :: myStatus
1363 +
1364 +  myStatus = 0
1365 +
1366 +  call complete_Shape_FF(myStatus)
1367 +
1368 +  status = myStatus
1369 +
1370 +  return
1371 + end subroutine completeShapeFF
1372 +

Diff Legend

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