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 1717 by chrisfen, Fri Nov 5 21:04:33 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
350      real (kind=dp), dimension(3), intent(inout) :: fpair
351 <    real (kind=dp) :: pot, vpair, sw
351 >    real (kind=dp) :: pot, vpair, sw, dswdr
352      real (kind=dp), dimension(9,nLocal) :: A
353      real (kind=dp), dimension(3,nLocal) :: f
354      real (kind=dp), dimension(3,nLocal) :: t
# Line 409 | Line 359 | contains  
359      integer :: l, m, lm, id1, id2, localError, function_type
360      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
361      real (kind=dp) :: coeff
362 +    real (kind=dp) :: pot_temp
363  
364      real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
365      real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
# Line 427 | Line 378 | contains  
378  
379      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
380  
381 +    real (kind=dp) :: sti2, stj2
382 +
383      real (kind=dp) :: proji, proji3, projj, projj3
384      real (kind=dp) :: cti, ctj, cpi, cpj, spi, spj
385      real (kind=dp) :: Phunc, sigma, s, eps, rtdenom, rt
# Line 485 | 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(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
442 +    real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
443 +    real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
444 +    real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
445 +
446      if (.not.haveShapeMap) then
447         call handleError("calc_shape", "NO SHAPEMAP!!!!")
448         return      
# Line 492 | 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 514 | Line 472 | contains  
472   #endif
473  
474      ! use the atid to find the shape type (st) for each atom:
517 –
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 559 | 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
566 <      
567 <       proji = sqrt(xi2 + yi2)
568 <       proji3 = proji*proji*proji
569 <      
523 >       zi2 = zi*zi            
524         cti = zi / rij
525 +
526 +       if (cti .gt. 1.0_dp) cti = 1.0_dp
527 +       if (cti .lt. -1.0_dp) cti = -1.0_dp
528 +
529         dctidx = - zi * xi / r3
530         dctidy = - zi * yi / r3
531         dctidz = 1.0d0 / rij - zi2 / r3
532 <       dctidux =  yi / rij
533 <       dctiduy = -xi / rij
534 <       dctiduz = 0.0d0
532 >       dctidux = - (zi * xi2) / r3
533 >       dctiduy = - (zi * yi2) / r3
534 >       dctiduz = zi / rij - (zi2 * zi) / r3
535 >
536 >       ! this is an attempt to try to truncate the singularity when
537 >       ! sin(theta) is near 0.0:
538 >
539 >       sti2 = 1.0_dp - cti*cti
540 >       if (dabs(sti2) .lt. 1.0d-12) then
541 >          proji = sqrt(rij * 1.0d-12)
542 >          dcpidx = 1.0d0 / proji
543 >          dcpidy = 0.0d0
544 >          dcpidux = xi / proji
545 >          dcpiduy = 0.0d0
546 >          dspidx = 0.0d0
547 >          dspidy = 1.0d0 / proji
548 >          dspidux = 0.0d0
549 >          dspiduy = yi / proji
550 >       else
551 >          proji = sqrt(xi2 + yi2)
552 >          proji3 = proji*proji*proji
553 >          dcpidx = 1.0d0 / proji - xi2 / proji3
554 >          dcpidy = - xi * yi / proji3
555 >          dcpidux = xi / proji - (xi2 * xi) / proji3
556 >          dcpiduy = - (xi * yi2) / proji3
557 >          dspidx = - xi * yi / proji3
558 >          dspidy = 1.0d0 / proji - yi2 / proji3
559 >          dspidux = - (yi * xi2) / proji3
560 >          dspiduy = yi / proji - (yi2 * yi) / proji3
561 >       endif
562        
563         cpi = xi / proji
579 –       dcpidx = 1.0d0 / proji - xi2 / proji3
580 –       dcpidy = - xi * yi / proji3
564         dcpidz = 0.0d0
565 <       dcpidux = xi * yi * zi / proji3
583 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
584 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
565 >       dcpiduz = 0.0d0
566        
567         spi = yi / proji
587 –       dspidx = - xi * yi / proji3
588 –       dspidy = 1.0d0 / proji - yi2 / proji3
568         dspidz = 0.0d0
569 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
591 <       dspiduy = xi * yi * zi / proji3
592 <       dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
569 >       dspiduz = 0.0d0
570  
571 <       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
572 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
571 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
572 >            ShapeMap%Shapes(st1)%bigL, LMAX, &
573 >            plm_i, dlm_i)
574  
575 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
575 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
576              CHEBYSHEV_TN, tm_i, dtm_i)
577 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
577 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
578              CHEBYSHEV_UN, um_i, dum_i)
579        
580         sigma_i = 0.0d0
# Line 645 | Line 623 | contains  
623               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
624            endif
625  
626 <          sigma_i = sigma_i + plm_i(l,m)*Phunc
626 >          sigma_i = sigma_i + plm_i(m,l)*Phunc
627 >
628 >          dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
629 >               Phunc * dlm_i(m,l) * dctidx
630 >          dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
631 >               Phunc * dlm_i(m,l) * dctidy
632 >          dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
633 >               Phunc * dlm_i(m,l) * dctidz
634            
635 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
636 <               Phunc * dlm_i(l,m) * dctidx
637 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
638 <               Phunc * dlm_i(l,m) * dctidy
639 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
640 <               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
635 >          dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
636 >               Phunc * dlm_i(m,l) * dctidux
637 >          dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
638 >               Phunc * dlm_i(m,l) * dctiduy
639 >          dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
640 >               Phunc * dlm_i(m,l) * dctiduz
641  
642         end do
643  
# Line 687 | Line 665 | contains  
665               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
666            endif
667  
668 <          s_i = s_i + plm_i(l,m)*Phunc
668 >          s_i = s_i + plm_i(m,l)*Phunc
669            
670 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
671 <               Phunc * dlm_i(l,m) * dctidx
672 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
673 <               Phunc * dlm_i(l,m) * dctidy
674 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
675 <               Phunc * dlm_i(l,m) * dctidz
670 >          dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
671 >               Phunc * dlm_i(m,l) * dctidx
672 >          dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
673 >               Phunc * dlm_i(m,l) * dctidy
674 >          dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
675 >               Phunc * dlm_i(m,l) * dctidz
676            
677 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
678 <               Phunc * dlm_i(l,m) * dctidux
679 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
680 <               Phunc * dlm_i(l,m) * dctiduy
681 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
682 <               Phunc * dlm_i(l,m) * dctiduz      
677 >          dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
678 >               Phunc * dlm_i(m,l) * dctidux
679 >          dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
680 >               Phunc * dlm_i(m,l) * dctiduy
681 >          dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
682 >               Phunc * dlm_i(m,l) * dctiduz      
683  
684         end do
685                
# Line 729 | Line 707 | contains  
707               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
708            endif
709  
710 <          eps_i = eps_i + plm_i(l,m)*Phunc
710 >          eps_i = eps_i + plm_i(m,l)*Phunc
711            
712 <          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
713 <               Phunc * dlm_i(l,m) * dctidx
714 <          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
715 <               Phunc * dlm_i(l,m) * dctidy
716 <          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
717 <               Phunc * dlm_i(l,m) * dctidz
712 >          depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
713 >               Phunc * dlm_i(m,l) * dctidx
714 >          depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
715 >               Phunc * dlm_i(m,l) * dctidy
716 >          depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
717 >               Phunc * dlm_i(m,l) * dctidz
718            
719 <          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
720 <               Phunc * dlm_i(l,m) * dctidux
721 <          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
722 <               Phunc * dlm_i(l,m) * dctiduy
723 <          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
724 <               Phunc * dlm_i(l,m) * dctiduz      
719 >          depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
720 >               Phunc * dlm_i(m,l) * dctidux
721 >          depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
722 >               Phunc * dlm_i(m,l) * dctiduy
723 >          depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
724 >               Phunc * dlm_i(m,l) * dctiduz      
725  
726         end do
727  
# Line 796 | Line 774 | contains  
774         xj2 = xj*xj
775         yj2 = yj*yj
776         zj2 = zj*zj
799 –      
800 –       projj = sqrt(xj2 + yj2)
801 –       projj3 = projj*projj*projj
802 –      
777         ctj = zj / rij
778 +      
779 +       if (ctj .gt. 1.0_dp) ctj = 1.0_dp
780 +       if (ctj .lt. -1.0_dp) ctj = -1.0_dp
781 +
782         dctjdx = - zj * xj / r3
783         dctjdy = - zj * yj / r3
784         dctjdz = 1.0d0 / rij - zj2 / r3
785 <       dctjdux =  yj / rij
786 <       dctjduy = -xj / rij
787 <       dctjduz = 0.0d0
785 >       dctjdux = - (zi * xj2) / r3
786 >       dctjduy = - (zj * yj2) / r3
787 >       dctjduz = zj / rij - (zj2 * zj) / r3
788        
789 <       cpj = xj / projj
790 <       dcpjdx = 1.0d0 / projj - xj2 / projj3
791 <       dcpjdy = - xj * yj / projj3
789 >       ! this is an attempt to try to truncate the singularity when
790 >       ! sin(theta) is near 0.0:
791 >
792 >       stj2 = 1.0_dp - ctj*ctj
793 >       if (dabs(stj2) .lt. 1.0d-12) then
794 >          projj = sqrt(rij * 1.0d-12)
795 >          dcpjdx = 1.0d0 / projj
796 >          dcpjdy = 0.0d0
797 >          dcpjdux = xj / projj
798 >          dcpjduy = 0.0d0
799 >          dspjdx = 0.0d0
800 >          dspjdy = 1.0d0 / projj
801 >          dspjdux = 0.0d0
802 >          dspjduy = yj / projj
803 >       else
804 >          projj = sqrt(xj2 + yj2)
805 >          projj3 = projj*projj*projj
806 >          dcpjdx = 1.0d0 / projj - xj2 / projj3
807 >          dcpjdy = - xj * yj / projj3
808 >          dcpjdux = xj / projj - (xj2 * xj) / projj3
809 >          dcpjduy = - (xj * yj2) / projj3
810 >          dspjdx = - xj * yj / projj3
811 >          dspjdy = 1.0d0 / projj - yj2 / projj3
812 >          dspjdux = - (yj * xj2) / projj3
813 >          dspjduy = yj / projj - (yj2 * yj) / projj3
814 >       endif
815 >
816 >       cpj = xj / projj
817         dcpjdz = 0.0d0
818 <       dcpjdux = xj * yj * zj / projj3
816 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
817 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
818 >       dcpjduz = 0.0d0
819        
820         spj = yj / projj
820 –       dspjdx = - xj * yj / projj3
821 –       dspjdy = 1.0d0 / projj - yj2 / projj3
821         dspjdz = 0.0d0
822 <       dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
823 <       dspjduy = xj * yj * zj / projj3
824 <       dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
822 >       dspjduz = 0.0d0
823 >
824 >
825 >       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
826 >       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
827 >       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
828 >            ShapeMap%Shapes(st2)%bigL, LMAX, &
829 >            plm_j, dlm_j)
830        
831 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
828 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
829 <      
830 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
831 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
832              CHEBYSHEV_TN, tm_j, dtm_j)
833 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
833 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
834              CHEBYSHEV_UN, um_j, dum_j)
835        
836         sigma_j = 0.0d0
# Line 878 | Line 879 | contains  
879               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
880            endif
881  
882 <          sigma_j = sigma_j + plm_j(l,m)*Phunc
882 >          sigma_j = sigma_j + plm_j(m,l)*Phunc
883            
884 <          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
885 <               Phunc * dlm_j(l,m) * dctjdx
886 <          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
887 <               Phunc * dlm_j(l,m) * dctjdy
888 <          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
889 <               Phunc * dlm_j(l,m) * dctjdz
884 >          dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
885 >               Phunc * dlm_j(m,l) * dctjdx
886 >          dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
887 >               Phunc * dlm_j(m,l) * dctjdy
888 >          dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
889 >               Phunc * dlm_j(m,l) * dctjdz
890            
891 <          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
892 <               Phunc * dlm_j(l,m) * dctjdux
893 <          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
894 <               Phunc * dlm_j(l,m) * dctjduy
895 <          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
896 <               Phunc * dlm_j(l,m) * dctjduz
891 >          dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
892 >               Phunc * dlm_j(m,l) * dctjdux
893 >          dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
894 >               Phunc * dlm_j(m,l) * dctjduy
895 >          dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
896 >               Phunc * dlm_j(m,l) * dctjduz
897  
898         end do
899  
# Line 920 | Line 921 | contains  
921               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
922            endif
923  
924 <          s_j = s_j + plm_j(l,m)*Phunc
924 >          s_j = s_j + plm_j(m,l)*Phunc
925            
926 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
927 <               Phunc * dlm_j(l,m) * dctjdx
928 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
929 <               Phunc * dlm_j(l,m) * dctjdy
930 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
931 <               Phunc * dlm_j(l,m) * dctjdz
926 >          dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
927 >               Phunc * dlm_j(m,l) * dctjdx
928 >          dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
929 >               Phunc * dlm_j(m,l) * dctjdy
930 >          dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
931 >               Phunc * dlm_j(m,l) * dctjdz
932            
933 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
934 <               Phunc * dlm_j(l,m) * dctjdux
935 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
936 <               Phunc * dlm_j(l,m) * dctjduy
937 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
938 <               Phunc * dlm_j(l,m) * dctjduz
933 >          dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
934 >               Phunc * dlm_j(m,l) * dctjdux
935 >          dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
936 >               Phunc * dlm_j(m,l) * dctjduy
937 >          dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
938 >               Phunc * dlm_j(m,l) * dctjduz
939  
940         end do
941  
# Line 962 | Line 963 | contains  
963               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
964            endif
965  
966 <          eps_j = eps_j + plm_j(l,m)*Phunc
966 >          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
967 >
968 >          eps_j = eps_j + plm_j(m,l)*Phunc
969            
970 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
971 <               Phunc * dlm_j(l,m) * dctjdx
972 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
973 <               Phunc * dlm_j(l,m) * dctjdy
974 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
975 <               Phunc * dlm_j(l,m) * dctjdz
970 >          depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
971 >               Phunc * dlm_j(m,l) * dctjdx
972 >          depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
973 >               Phunc * dlm_j(m,l) * dctjdy
974 >          depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
975 >               Phunc * dlm_j(m,l) * dctjdz
976            
977 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
978 <               Phunc * dlm_j(l,m) * dctjdux
979 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
980 <               Phunc * dlm_j(l,m) * dctjduy
981 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
982 <               Phunc * dlm_j(l,m) * dctjduz
977 >          depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
978 >               Phunc * dlm_j(m,l) * dctjdux
979 >          depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
980 >               Phunc * dlm_j(m,l) * dctjduy
981 >          depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
982 >               Phunc * dlm_j(m,l) * dctjduz
983  
984         end do
985  
# Line 1032 | Line 1035 | contains  
1035      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1036      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1037      
1038 + !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1039 + !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1040 + !!$
1041 + !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1042 + !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1043 + !!$
1044 + !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1045 +
1046      rtdenom = rij-sigma+s
1047      rt = s / rtdenom
1048  
# Line 1056 | Line 1067 | contains  
1067      rt12 = rt6*rt6
1068      rt126 = rt12 - rt6
1069  
1070 +    pot_temp = 4.0d0 * eps * rt126
1071 +
1072 +    vpair = vpair + pot_temp
1073      if (do_pot) then
1074   #ifdef IS_MPI
1075 <       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1076 <       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1075 >       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1076 >       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1077   #else
1078 <       pot = pot + 4.0d0*eps*rt126*sw
1078 >       pot = pot + pot_temp*sw
1079   #endif
1080      endif
1081 +
1082 + !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1083      
1084      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1085      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
# Line 1082 | Line 1098 | contains  
1098      ! do the torques first since they are easy:
1099      ! remember that these are still in the body fixed axes
1100  
1085 –    txi = dvduxi * sw
1086 –    tyi = dvduyi * sw
1087 –    tzi = dvduzi * sw
1101  
1102 <    txj = dvduxj * sw
1103 <    tyj = dvduyj * sw
1104 <    tzj = dvduzj * sw
1102 > !!$    write(*,*) 'sw = ', sw
1103 > !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1104 > !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1105 > !!$
1106 >    txi =  (dvduzi - dvduyi) * sw
1107 >    tyi =  (dvduxi - dvduzi) * sw
1108 >    tzi =  (dvduyi - dvduxi) * sw
1109  
1110 +    txj = (dvduzj - dvduyj) * sw
1111 +    tyj = (dvduxj - dvduzj) * sw
1112 +    tzj = (dvduyj - dvduxj) * sw
1113 +
1114 + !!$    txi = -dvduxi * sw
1115 + !!$    tyi = -dvduyi * sw
1116 + !!$    tzi = -dvduzi * sw
1117 + !!$
1118 + !!$    txj = dvduxj * sw
1119 + !!$    tyj = dvduyj * sw
1120 + !!$    tzj = dvduzj * sw
1121 +
1122 +    write(*,*) 't1 = ', txi, tyi, tzi
1123 +    write(*,*) 't2 = ', txj, tyj, tzj
1124 +
1125      ! go back to lab frame using transpose of rotation matrix:
1126      
1127   #ifdef IS_MPI
# Line 1153 | Line 1185 | contains  
1185      fyji = -fyjj
1186      fzji = -fzjj
1187  
1188 <    fxradial = fxii + fxji
1189 <    fyradial = fyii + fyji
1190 <    fzradial = fzii + fzji
1188 >    fxradial = 0.5_dp * (fxii + fxji)
1189 >    fyradial = 0.5_dp * (fyii + fyji)
1190 >    fzradial = 0.5_dp * (fzii + fzji)
1191  
1192   #ifdef IS_MPI
1193      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
# Line 1190 | Line 1222 | contains  
1222         fpair(3) = fpair(3) + fzradial
1223        
1224      endif
1225 <    
1225 >
1226    end subroutine do_shape_pair
1227      
1228 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1229 <    
1228 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1229 >
1230      ! Purpose: Compute the associated Legendre functions
1231      !          Plm(x) and their derivatives Plm'(x)
1232      ! Input :  x  --- Argument of Plm(x)
# Line 1211 | Line 1243 | contains  
1243      ! The original Fortran77 codes can be found here:
1244      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1245      
1246 <    real (kind=8), intent(in) :: x
1246 >    real (kind=dp), intent(in) :: x
1247      integer, intent(in) :: l, m, lmax
1248 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1248 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1249      integer :: i, j, ls
1250 <    real (kind=8) :: xq, xs
1250 >    real (kind=dp) :: xq, xs
1251  
1252      ! zero out both arrays:
1253      DO I = 0, m
1254         DO J = 0, l
1255 <          PLM(J,I) = 0.0D0
1256 <          DLM(J,I) = 0.0D0
1255 >          PLM(J,I) = 0.0_dp
1256 >          DLM(J,I) = 0.0_dp
1257         end DO
1258      end DO
1259  
# Line 1254 | Line 1286 | contains  
1286      DO I = 1, l
1287         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1288      enddo
1289 <    
1289 >
1290      DO I = 0, l
1291         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1292      enddo
1293 <    
1293 >
1294      DO I = 0, l
1295         DO J = I+2, m
1296            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1297                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1298         end DO
1299      end DO
1300 <  
1300 >
1301      DLM(0, 0)=0.0D0
1270 –    
1302      DO J = 1, m
1303         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1304      end DO
1305 <    
1305 >
1306      DO I = 1, l
1307         DO J = I, m
1308            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1309         end DO
1310      end DO
1311 <    
1311 >
1312      RETURN
1313    END SUBROUTINE Associated_Legendre
1314  
1315  
1316 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1316 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1317    
1318      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1319      !          or Ln(x) or Hn(x), and their derivatives
# Line 1304 | Line 1335 | contains  
1335      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1336    
1337      real(kind=8), intent(in) :: x
1338 <    integer, intent(in):: m
1338 >    integer, intent(in):: m, mmax
1339      integer, intent(in):: function_type
1340 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1340 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1341    
1342      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1343      integer :: k
# Line 1350 | Line 1381 | contains  
1381         DY0 = DY1
1382         DY1 = DYN
1383      end DO
1384 +
1385 +
1386      RETURN
1387      
1388    end subroutine Orthogonal_Polynomial
# Line 1361 | Line 1394 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1394       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1395       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1396       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1397 <     myAtid, status)
1397 >     myATID, status)
1398  
1399    use definitions
1400    use shapes, only: newShapeType
# Line 1370 | Line 1403 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1403    integer :: nRangeFuncs
1404    integer :: nStrengthFuncs
1405    integer :: status
1406 <  integer :: myAtid
1406 >  integer :: myATID
1407    
1408    integer, dimension(nContactFuncs) :: ContactFuncLValue          
1409    integer, dimension(nContactFuncs) :: ContactFuncMValue          
# Line 1390 | Line 1423 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1423         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1424         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1425         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1426 <       myAtid, status)
1426 >       myATID, status)
1427  
1428    return
1429   end subroutine makeShape
1430 +
1431 + subroutine completeShapeFF(status)
1432 +
1433 +  use shapes, only: complete_Shape_FF
1434 +
1435 +  integer, intent(out)  :: status
1436 +  integer :: myStatus
1437 +
1438 +  myStatus = 0
1439 +
1440 +  call complete_Shape_FF(myStatus)
1441 +
1442 +  status = myStatus
1443 +
1444 +  return
1445 + end subroutine completeShapeFF
1446 +

Diff Legend

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