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 1647 by chrisfen, Tue Oct 26 18:03:12 2004 UTC vs.
Revision 1707 by gezelter, Thu Nov 4 16:20:37 2004 UTC

# Line 62 | Line 62 | module shapes
62    type(ShapeList), save :: ShapeMap
63  
64    integer :: lmax
65 –  real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
66 –  real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
67 –  real (kind=dp), allocatable, dimension(:) :: tm_j, dtm_j, um_j, dum_j
65  
66   contains  
67  
# Line 73 | 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 115 | 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 128 | 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 188 | 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 253 | 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 287 | Line 286 | contains  
286         return
287      endif
288  
289 +    return
290 +
291    end subroutine allocateShape
292      
293    subroutine complete_Shape_FF(status)
# Line 310 | 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_LennardJones", thisProperty)
318        
# Line 339 | Line 341 | contains  
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 372 | Line 377 | contains  
377  
378      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
379  
380 +    real (kind=dp) :: sti2, stj2
381 +
382      real (kind=dp) :: proji, proji3, projj, projj3
383      real (kind=dp) :: cti, ctj, cpi, cpj, spi, spj
384      real (kind=dp) :: Phunc, sigma, s, eps, rtdenom, rt
# Line 430 | Line 437 | contains  
437      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
438      real (kind=dp) :: fxradial, fyradial, fzradial
439  
440 +    real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
441 +    real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
442 +    real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
443 +    real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
444 +
445      if (.not.haveShapeMap) then
446         call handleError("calc_shape", "NO SHAPEMAP!!!!")
447         return      
448      endif
449 +
450 +    write(*,*) rij, r2, d(1), d(2), d(3)
451 +    write(*,*) 'before, atom1, 2 = ', atom1, atom2
452 +    write(*,*) 'f1 = ', f(1,atom1), f(2,atom1), f(3,atom1)
453 +    write(*,*) 'f2 = ', f(1,atom2), f(2,atom2), f(3,atom2)
454 +    write(*,*) 't1 = ', t(1,atom1), t(2,atom1), t(3,atom1)
455 +    write(*,*) 't2 = ', t(1,atom2), t(2,atom2), t(3,atom2)
456      
457      !! We assume that the rotation matrices have already been calculated
458      !! and placed in the A array.
459 <        
459 >
460      r3 = r2*rij
461      r5 = r3*r2
462      
# Line 459 | Line 478 | contains  
478   #endif
479  
480      ! use the atid to find the shape type (st) for each atom:
462 –
481      st1 = ShapeMap%atidToShape(atid1)
482      st2 = ShapeMap%atidToShape(atid2)
483 <    
483 >
484      if (ShapeMap%Shapes(st1)%isLJ) then
485 +
486         sigma_i = ShapeMap%Shapes(st1)%sigma
487         s_i = ShapeMap%Shapes(st1)%sigma
488         eps_i = ShapeMap%Shapes(st1)%epsilon
# Line 504 | Line 523 | contains  
523         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
524        
525   #endif
526 <      
526 >
527         xi2 = xi*xi
528         yi2 = yi*yi
529 <       zi2 = zi*zi
511 <      
512 <       proji = sqrt(xi2 + yi2)
513 <       proji3 = proji*proji*proji
514 <      
529 >       zi2 = zi*zi            
530         cti = zi / rij
531 +
532 +       if (cti .gt. 1.0_dp) cti = 1.0_dp
533 +       if (cti .lt. -1.0_dp) cti = -1.0_dp
534 +
535         dctidx = - zi * xi / r3
536         dctidy = - zi * yi / r3
537         dctidz = 1.0d0 / rij - zi2 / r3
538 <       dctidux =  yi / rij
539 <       dctiduy = -xi / rij
538 >       dctidux =  yi / rij + (zi * yi) / r3
539 >       dctiduy = -xi / rij - (zi * xi) / r3
540         dctiduz = 0.0d0
541 +
542 +       ! this is an attempt to try to truncate the singularity when
543 +       ! sin(theta) is near 0.0:
544 +
545 +       sti2 = 1.0_dp - cti*cti
546 +       if (dabs(sti2) .lt. 1.0d-12) then
547 +          proji = sqrt(rij * 1.0d-12)
548 +          dcpidx = 1.0d0 / proji
549 +          dcpidy = 0.0d0
550 +          dcpidux = 0.0d0
551 +          dcpiduy = zi / proji
552 +          dspidx = 0.0d0
553 +          dspidy = 1.0d0 / proji
554 +          dspidux = -zi / proji
555 +          dspiduy = 0.0d0
556 +       else
557 +          proji = sqrt(xi2 + yi2)
558 +          proji3 = proji*proji*proji
559 +          dcpidx = 1.0d0 / proji - xi2 / proji3
560 +          dcpidy = - xi * yi / proji3
561 +          dcpidux = xi * yi * zi / proji3
562 +          dcpiduy = zi / proji - xi2 * zi / proji3
563 +          dspidx = - xi * yi / proji3
564 +          dspidy = 1.0d0 / proji - yi2 / proji3
565 +          dspidux = -zi / proji + yi2 * zi / proji3
566 +          dspiduy = - xi * yi * zi / proji3
567 +       endif
568        
569         cpi = xi / proji
524 –       dcpidx = 1.0d0 / proji - xi2 / proji3
525 –       dcpidy = - xi * yi / proji3
570         dcpidz = 0.0d0
571 <       dcpidux = xi * yi * zi / proji3
528 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
529 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
571 >       dcpiduz = -yi / proji
572        
573         spi = yi / proji
532 –       dspidx = - xi * yi / proji3
533 –       dspidy = 1.0d0 / proji - yi2 / proji3
574         dspidz = 0.0d0
575 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
576 <       dspiduy = xi * yi * zi / proji3
537 <       dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
575 >       dspiduz = xi / proji
576 >       write(*,*) 'before lmloop', cpi, dcpidx, dcpidux
577  
578 <       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
579 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
578 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
579 >            ShapeMap%Shapes(st1)%bigL, LMAX, &
580 >            plm_i, dlm_i)
581  
582 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
582 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
583              CHEBYSHEV_TN, tm_i, dtm_i)
584 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
584 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
585              CHEBYSHEV_UN, um_i, dum_i)
586        
587         sigma_i = 0.0d0
# Line 590 | Line 630 | contains  
630               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
631            endif
632  
633 <          sigma_i = sigma_i + plm_i(l,m)*Phunc
633 >          sigma_i = sigma_i + plm_i(m,l)*Phunc
634 >
635 >          dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
636 >               Phunc * dlm_i(m,l) * dctidx
637 >          dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
638 >               Phunc * dlm_i(m,l) * dctidy
639 >          dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
640 >               Phunc * dlm_i(m,l) * dctidz
641            
642 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
643 <               Phunc * dlm_i(l,m) * dctidx
644 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
645 <               Phunc * dlm_i(l,m) * dctidy
646 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
647 <               Phunc * dlm_i(l,m) * dctidz
601 <          
602 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
603 <               Phunc * dlm_i(l,m) * dctidux
604 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
605 <               Phunc * dlm_i(l,m) * dctiduy
606 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
607 <               Phunc * dlm_i(l,m) * dctiduz
642 >          dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
643 >               Phunc * dlm_i(m,l) * dctidux
644 >          dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
645 >               Phunc * dlm_i(m,l) * dctiduy
646 >          dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
647 >               Phunc * dlm_i(m,l) * dctiduz
648  
649         end do
650  
# Line 614 | Line 654 | contains  
654            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
655            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
656            
657 +         write(*,*) 'in lm loop a', coeff, dtm_i(m), dcpidx
658 +
659 +
660            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
661               Phunc = coeff * tm_i(m)
662               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 632 | Line 675 | contains  
675               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
676            endif
677  
678 <          s_i = s_i + plm_i(l,m)*Phunc
678 >          s_i = s_i + plm_i(m,l)*Phunc
679 >          
680 >
681 >          write(*,*) 'in lm loop ', dsidx, plm_i(m,l), dPhuncdX, Phunc, dlm_i(m,l), dctidx
682 >          dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
683 >               Phunc * dlm_i(m,l) * dctidx
684 >          dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
685 >               Phunc * dlm_i(m,l) * dctidy
686 >          dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
687 >               Phunc * dlm_i(m,l) * dctidz
688            
689 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
690 <               Phunc * dlm_i(l,m) * dctidx
691 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
692 <               Phunc * dlm_i(l,m) * dctidy
693 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
694 <               Phunc * dlm_i(l,m) * dctidz
643 <          
644 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
645 <               Phunc * dlm_i(l,m) * dctidux
646 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
647 <               Phunc * dlm_i(l,m) * dctiduy
648 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
649 <               Phunc * dlm_i(l,m) * dctiduz      
689 >          dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
690 >               Phunc * dlm_i(m,l) * dctidux
691 >          dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
692 >               Phunc * dlm_i(m,l) * dctiduy
693 >          dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
694 >               Phunc * dlm_i(m,l) * dctiduz      
695  
696         end do
697                
# Line 674 | Line 719 | contains  
719               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
720            endif
721  
722 <          eps_i = eps_i + plm_i(l,m)*Phunc
722 >          eps_i = eps_i + plm_i(m,l)*Phunc
723            
724 <          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
725 <               Phunc * dlm_i(l,m) * dctidx
726 <          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
727 <               Phunc * dlm_i(l,m) * dctidy
728 <          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
729 <               Phunc * dlm_i(l,m) * dctidz
724 >          depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
725 >               Phunc * dlm_i(m,l) * dctidx
726 >          depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
727 >               Phunc * dlm_i(m,l) * dctidy
728 >          depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
729 >               Phunc * dlm_i(m,l) * dctidz
730            
731 <          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
732 <               Phunc * dlm_i(l,m) * dctidux
733 <          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
734 <               Phunc * dlm_i(l,m) * dctiduy
735 <          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
736 <               Phunc * dlm_i(l,m) * dctiduz      
731 >          depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
732 >               Phunc * dlm_i(m,l) * dctidux
733 >          depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
734 >               Phunc * dlm_i(m,l) * dctiduy
735 >          depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
736 >               Phunc * dlm_i(m,l) * dctiduz      
737  
738         end do
739  
# Line 741 | Line 786 | contains  
786         xj2 = xj*xj
787         yj2 = yj*yj
788         zj2 = zj*zj
744 –      
745 –       projj = sqrt(xj2 + yj2)
746 –       projj3 = projj*projj*projj
747 –      
789         ctj = zj / rij
790 +      
791 +       if (ctj .gt. 1.0_dp) ctj = 1.0_dp
792 +       if (ctj .lt. -1.0_dp) ctj = -1.0_dp
793 +
794         dctjdx = - zj * xj / r3
795         dctjdy = - zj * yj / r3
796         dctjdz = 1.0d0 / rij - zj2 / r3
797 <       dctjdux =  yj / rij
798 <       dctjduy = -xj / rij
797 >       dctjdux =  yj / rij + (zj * yj) / r3
798 >       dctjduy = -xj / rij - (zj * xj) / r3
799         dctjduz = 0.0d0
800        
801 <       cpj = xj / projj
802 <       dcpjdx = 1.0d0 / projj - xj2 / projj3
803 <       dcpjdy = - xj * yj / projj3
804 <       dcpjdz = 0.0d0
805 <       dcpjdux = xj * yj * zj / projj3
806 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
807 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
801 >       ! this is an attempt to try to truncate the singularity when
802 >       ! sin(theta) is near 0.0:
803 >
804 >       stj2 = 1.0_dp - ctj*ctj
805 >       if (dabs(stj2) .lt. 1.0d-12) then
806 >          projj = sqrt(rij * 1.0d-12)
807 >          dcpjdx = 1.0d0 / projj
808 >          dcpjdy = 0.0d0
809 >          dcpjdux = 0.0d0
810 >          dcpjduy = zj / projj
811 >          dspjdx = 0.0d0
812 >          dspjdy = 1.0d0 / projj
813 >          dspjdux = -zj / projj
814 >          dspjduy = 0.0d0
815 >       else
816 >          projj = sqrt(xj2 + yj2)
817 >          projj3 = projj*projj*projj
818 >          dcpjdx = 1.0d0 / projj - xj2 / projj3
819 >          dcpjdy = - xj * yj / projj3
820 >          dcpjdux = xj * yj * zj / projj3
821 >          dcpjduy = zj / projj - xj2 * zj / projj3
822 >          dspjdx = - xj * yj / projj3
823 >          dspjdy = 1.0d0 / projj - yj2 / projj3
824 >          dspjdux = -zj / projj + yj2 * zj / projj3
825 >          dspjduy = - xj * yj * zj / projj3
826 >       endif
827 >
828 >       cpj = xj / projj
829 >       dcpjdz = 0.0d0
830 >       dcpjduz = -yj / projj
831        
832         spj = yj / projj
765 –       dspjdx = - xj * yj / projj3
766 –       dspjdy = 1.0d0 / projj - yj2 / projj3
833         dspjdz = 0.0d0
834 <       dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
835 <       dspjduy = xj * yj * zj / projj3
836 <       dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
834 >       dspjduz = xj / projj
835 >
836 >       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
837 >            ShapeMap%Shapes(st2)%bigL, LMAX, &
838 >            plm_j, dlm_j)
839        
840 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
773 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
774 <      
775 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
840 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
841              CHEBYSHEV_TN, tm_j, dtm_j)
842 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
842 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
843              CHEBYSHEV_UN, um_j, dum_j)
844        
845         sigma_j = 0.0d0
# Line 823 | Line 888 | contains  
888               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
889            endif
890  
891 <          sigma_j = sigma_j + plm_j(l,m)*Phunc
891 >          sigma_j = sigma_j + plm_j(m,l)*Phunc
892            
893 <          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
894 <               Phunc * dlm_j(l,m) * dctjdx
895 <          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
896 <               Phunc * dlm_j(l,m) * dctjdy
897 <          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
898 <               Phunc * dlm_j(l,m) * dctjdz
893 >          dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
894 >               Phunc * dlm_j(m,l) * dctjdx
895 >          dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
896 >               Phunc * dlm_j(m,l) * dctjdy
897 >          dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
898 >               Phunc * dlm_j(m,l) * dctjdz
899            
900 <          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
901 <               Phunc * dlm_j(l,m) * dctjdux
902 <          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
903 <               Phunc * dlm_j(l,m) * dctjduy
904 <          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
905 <               Phunc * dlm_j(l,m) * dctjduz
900 >          dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
901 >               Phunc * dlm_j(m,l) * dctjdux
902 >          dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
903 >               Phunc * dlm_j(m,l) * dctjduy
904 >          dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
905 >               Phunc * dlm_j(m,l) * dctjduz
906  
907         end do
908  
# Line 865 | Line 930 | contains  
930               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
931            endif
932  
933 <          s_j = s_j + plm_j(l,m)*Phunc
933 >          s_j = s_j + plm_j(m,l)*Phunc
934            
935 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
936 <               Phunc * dlm_j(l,m) * dctjdx
937 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
938 <               Phunc * dlm_j(l,m) * dctjdy
939 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
940 <               Phunc * dlm_j(l,m) * dctjdz
935 >          dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
936 >               Phunc * dlm_j(m,l) * dctjdx
937 >          dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
938 >               Phunc * dlm_j(m,l) * dctjdy
939 >          dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
940 >               Phunc * dlm_j(m,l) * dctjdz
941            
942 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
943 <               Phunc * dlm_j(l,m) * dctjdux
944 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
945 <               Phunc * dlm_j(l,m) * dctjduy
946 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
947 <               Phunc * dlm_j(l,m) * dctjduz
942 >          dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
943 >               Phunc * dlm_j(m,l) * dctjdux
944 >          dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
945 >               Phunc * dlm_j(m,l) * dctjduy
946 >          dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
947 >               Phunc * dlm_j(m,l) * dctjduz
948  
949         end do
950  
# Line 907 | Line 972 | contains  
972               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
973            endif
974  
975 <          eps_j = eps_j + plm_j(l,m)*Phunc
975 >          eps_j = eps_j + plm_j(m,l)*Phunc
976            
977 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
978 <               Phunc * dlm_j(l,m) * dctjdx
979 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
980 <               Phunc * dlm_j(l,m) * dctjdy
981 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
982 <               Phunc * dlm_j(l,m) * dctjdz
977 >          depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
978 >               Phunc * dlm_j(m,l) * dctjdx
979 >          depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
980 >               Phunc * dlm_j(m,l) * dctjdy
981 >          depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
982 >               Phunc * dlm_j(m,l) * dctjdz
983            
984 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
985 <               Phunc * dlm_j(l,m) * dctjdux
986 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
987 <               Phunc * dlm_j(l,m) * dctjduy
988 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
989 <               Phunc * dlm_j(l,m) * dctjduz
984 >          depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
985 >               Phunc * dlm_j(m,l) * dctjdux
986 >          depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
987 >               Phunc * dlm_j(m,l) * dctjduy
988 >          depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
989 >               Phunc * dlm_j(m,l) * dctjduz
990  
991         end do
992  
# Line 963 | Line 1028 | contains  
1028  
1029      eps = sqrt(eps_i * eps_j)
1030  
1031 +    write(*,*) 'sigma, s, eps = ', sigma, s, eps
1032 +
1033      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1034      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1035      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 978 | Line 1045 | contains  
1045      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1046      
1047      rtdenom = rij-sigma+s
1048 +
1049 +    write(*,*) 'rtdenom = ', rtdenom, ' sw = ', sw
1050      rt = s / rtdenom
1051  
1052 +    write(*,*) 'john' , dsdxi, rt, drdxi, dsigmadxi, rtdenom
1053 +    write(*,*) 'bigboot', dsduzj, rt, drduzj, dsigmaduzj, rtdenom
1054 +
1055 +
1056      drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1057      drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1058      drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
# Line 1010 | Line 1083 | contains  
1083   #endif
1084      endif
1085      
1086 +    write(*,*) 'drtdxi = ', drtdxi, drtdyi
1087 +    write(*,*) 'depsdxi = ', depsdxi, depsdyi
1088 +
1089      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1090      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1091      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1017 | Line 1093 | contains  
1093      dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1094      dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1095  
1096 +    write(*,*) 'drtduzj = ', drtduzj, depsduzj
1097 +
1098      dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1099      dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1100      dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
# Line 1027 | Line 1105 | contains  
1105      ! do the torques first since they are easy:
1106      ! remember that these are still in the body fixed axes
1107  
1108 +    write(*,*) 'dvdx = ', dvdxi, dvdyi, dvdzi
1109 +    write(*,*) 'dvdx = ', dvdxj, dvdyj, dvdzj
1110 +    write(*,*) 'dvdu = ', dvduxi, dvduyi, dvduzi
1111 +    write(*,*) 'dvdu = ', dvduxj, dvduyj, dvduzj
1112 +
1113      txi = dvduxi * sw
1114      tyi = dvduyi * sw
1115      tzi = dvduzi * sw
# Line 1136 | Line 1219 | contains  
1219        
1220      endif
1221      
1222 +    write(*,*) 'f1 = ', f(1,atom1), f(2,atom1), f(3,atom1)
1223 +    write(*,*) 'f2 = ', f(1,atom2), f(2,atom2), f(3,atom2)
1224 +    write(*,*) 't1 = ', t(1,atom1), t(2,atom1), t(3,atom1)
1225 +    write(*,*) 't2 = ', t(1,atom2), t(2,atom2), t(3,atom2)
1226 +
1227 +
1228    end subroutine do_shape_pair
1229      
1230 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1231 <    
1230 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1231 >
1232      ! Purpose: Compute the associated Legendre functions
1233      !          Plm(x) and their derivatives Plm'(x)
1234      ! Input :  x  --- Argument of Plm(x)
# Line 1156 | Line 1245 | contains  
1245      ! The original Fortran77 codes can be found here:
1246      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1247      
1248 <    real (kind=8), intent(in) :: x
1248 >    real (kind=dp), intent(in) :: x
1249      integer, intent(in) :: l, m, lmax
1250 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1250 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1251      integer :: i, j, ls
1252 <    real (kind=8) :: xq, xs
1252 >    real (kind=dp) :: xq, xs
1253  
1254      ! zero out both arrays:
1255      DO I = 0, m
1256         DO J = 0, l
1257 <          PLM(J,I) = 0.0D0
1258 <          DLM(J,I) = 0.0D0
1257 >          PLM(J,I) = 0.0_dp
1258 >          DLM(J,I) = 0.0_dp
1259         end DO
1260      end DO
1261  
# Line 1199 | Line 1288 | contains  
1288      DO I = 1, l
1289         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1290      enddo
1291 <    
1291 >
1292      DO I = 0, l
1293         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1294      enddo
1295 <    
1295 >
1296      DO I = 0, l
1297         DO J = I+2, m
1298            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1299                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1300         end DO
1301      end DO
1302 <  
1302 >
1303      DLM(0, 0)=0.0D0
1215 –    
1304      DO J = 1, m
1305         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1306      end DO
1307 <    
1307 >
1308      DO I = 1, l
1309         DO J = I, m
1310            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1311         end DO
1312      end DO
1313 <    
1313 >
1314      RETURN
1315    END SUBROUTINE Associated_Legendre
1316  
1317  
1318 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1318 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1319    
1320      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1321      !          or Ln(x) or Hn(x), and their derivatives
# Line 1249 | Line 1337 | contains  
1337      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1338    
1339      real(kind=8), intent(in) :: x
1340 <    integer, intent(in):: m
1340 >    integer, intent(in):: m, mmax
1341      integer, intent(in):: function_type
1342 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1342 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1343    
1344      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1345      integer :: k
# Line 1295 | Line 1383 | contains  
1383         DY0 = DY1
1384         DY1 = DYN
1385      end DO
1386 +
1387 +
1388      RETURN
1389      
1390    end subroutine Orthogonal_Polynomial
# Line 1306 | Line 1396 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1396       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1397       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1398       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1399 <     myAtid, status)
1399 >     myATID, status)
1400  
1401    use definitions
1402    use shapes, only: newShapeType
# Line 1315 | Line 1405 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1405    integer :: nRangeFuncs
1406    integer :: nStrengthFuncs
1407    integer :: status
1408 <  integer :: myAtid
1408 >  integer :: myATID
1409    
1410    integer, dimension(nContactFuncs) :: ContactFuncLValue          
1411    integer, dimension(nContactFuncs) :: ContactFuncMValue          
# Line 1335 | Line 1425 | subroutine makeShape(nContactFuncs, ContactFuncLValue,
1425         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1426         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1427         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1428 <       myAtid, status)
1428 >       myATID, status)
1429  
1430    return
1431   end subroutine makeShape

Diff Legend

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