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 1689 by chrisfen, Fri Oct 29 22:28:45 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 380 | 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 438 | 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(LMAX,MMAX), dlm_i(LMAX,MMAX)
441 <    real (kind=dp) :: plm_j(LMAX,MMAX), dlm_j(LMAX,MMAX)
442 <    real (kind=dp) :: tm_i(MMAX), dtm_i(MMAX), um_i(MMAX), dum_i(MMAX)
443 <    real (kind=dp) :: tm_j(MMAX), dtm_j(MMAX), um_j(MMAX), dum_j(MMAX)
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.
# Line 520 | Line 526 | contains  
526  
527         xi2 = xi*xi
528         yi2 = yi*yi
529 <       zi2 = zi*zi
524 <      
525 <       proji = sqrt(xi2 + yi2)
526 <       proji3 = proji*proji*proji
527 <      
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
538 –       dcpidx = 1.0d0 / proji - xi2 / proji3
539 –       dcpidy = - xi * yi / proji3
570         dcpidz = 0.0d0
571 <       dcpidux = xi * yi * zi / proji3
542 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
543 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
571 >       dcpiduz = -yi / proji
572        
573         spi = yi / proji
546 –       dspidx = - xi * yi / proji3
547 –       dspidy = 1.0d0 / proji - yi2 / proji3
574         dspidz = 0.0d0
575 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
576 <       dspiduy = xi * yi * zi / proji3
551 <       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, &
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, MMAX, &
# Line 588 | Line 613 | contains  
613            function_type = ShapeMap%Shapes(st1)%ContactFunctionType(lm)
614  
615            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
591 –           !  write(*,*) tm_i(m), ' is tm_i'
616               Phunc = coeff * tm_i(m)
617               dPhuncdX = coeff * dtm_i(m) * dcpidx
618               dPhuncdY = coeff * dtm_i(m) * dcpidy
# Line 606 | 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
634 <          write(*,*) plm_i(l,m), l, m
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
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 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
643 <               Phunc * dlm_i(l,m) * dctidux
644 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
645 <               Phunc * dlm_i(l,m) * dctiduy
646 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
647 <               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 630 | 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 648 | 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
659 <          
660 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
661 <               Phunc * dlm_i(l,m) * dctidux
662 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
663 <               Phunc * dlm_i(l,m) * dctiduy
664 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
665 <               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 689 | Line 718 | contains  
718               dPhuncdUy = coeff*(spi * dum_i(m-1)*dcpiduy + dspiduy *um_i(m-1))
719               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
720            endif
721 <          !write(*,*) eps_i, plm_i(l,m), Phunc
722 <          eps_i = eps_i + plm_i(l,m)*Phunc
721 >
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 757 | Line 786 | contains  
786         xj2 = xj*xj
787         yj2 = yj*yj
788         zj2 = zj*zj
760 –      
761 –       projj = sqrt(xj2 + yj2)
762 –       projj3 = projj*projj*projj
763 –      
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 +       ! 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
773 –       dcpjdx = 1.0d0 / projj - xj2 / projj3
774 –       dcpjdy = - xj * yj / projj3
829         dcpjdz = 0.0d0
830 <       dcpjdux = xj * yj * zj / projj3
777 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
778 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
830 >       dcpjduz = -yj / projj
831        
832         spj = yj / projj
781 –       dspjdx = - xj * yj / projj3
782 –       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)
837 <      
788 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
789 <            ShapeMap%Shapes(st2)%bigM, LMAX, &
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 Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
# Line 840 | 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 882 | 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 924 | 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 977 | Line 1025 | contains  
1025      dsduxj = 0.5*dsjdux
1026      dsduyj = 0.5*dsjduy
1027      dsduzj = 0.5*dsjduz
1028 <    !write(*,*) eps_i, eps_j
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 995 | 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 1027 | 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 1034 | 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 1044 | 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 1153 | 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)        
# Line 1311 | Line 1383 | contains  
1383         DY0 = DY1
1384         DY1 = DYN
1385      end DO
1386 +
1387 +
1388      RETURN
1389      
1390    end subroutine Orthogonal_Polynomial

Diff Legend

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