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 1717 by chrisfen, Fri Nov 5 21:04:33 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 351 | Line 348 | contains  
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 362 | 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 379 | Line 377 | contains  
377      real (kind=dp) :: depsjdux, depsjduy, depsjduz
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
# Line 438 | Line 438 | contains  
438      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
439      real (kind=dp) :: fxradial, fyradial, fzradial
440  
441 <    real (kind=dp) :: plm_i(LMAX,MMAX), dlm_i(LMAX,MMAX)
442 <    real (kind=dp) :: plm_j(LMAX,MMAX), dlm_j(LMAX,MMAX)
443 <    real (kind=dp) :: tm_i(MMAX), dtm_i(MMAX), um_i(MMAX), dum_i(MMAX)
444 <    real (kind=dp) :: tm_j(MMAX), dtm_j(MMAX), um_j(MMAX), dum_j(MMAX)
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!!!!")
# Line 520 | Line 520 | contains  
520  
521         xi2 = xi*xi
522         yi2 = yi*yi
523 <       zi2 = zi*zi
524 <      
525 <       proji = sqrt(xi2 + yi2)
526 <       proji3 = proji*proji*proji
527 <      
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
538 –       dcpidx = 1.0d0 / proji - xi2 / proji3
539 –       dcpidy = - xi * yi / proji3
564         dcpidz = 0.0d0
565 <       dcpidux = xi * yi * zi / proji3
542 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
543 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
565 >       dcpiduz = 0.0d0
566        
567         spi = yi / proji
546 –       dspidx = - xi * yi / proji3
547 –       dspidy = 1.0d0 / proji - yi2 / proji3
568         dspidz = 0.0d0
569 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
550 <       dspiduy = xi * yi * zi / proji3
551 <       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, &
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, MMAX, &
# Line 588 | Line 606 | contains  
606            function_type = ShapeMap%Shapes(st1)%ContactFunctionType(lm)
607  
608            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
591 –           !  write(*,*) tm_i(m), ' is tm_i'
609               Phunc = coeff * tm_i(m)
610               dPhuncdX = coeff * dtm_i(m) * dcpidx
611               dPhuncdY = coeff * dtm_i(m) * dcpidy
# Line 606 | 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
627 <          write(*,*) plm_i(l,m), l, m
628 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
629 <               Phunc * dlm_i(l,m) * dctidx
630 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
631 <               Phunc * dlm_i(l,m) * dctidy
632 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
633 <               Phunc * dlm_i(l,m) * dctidz
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 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
636 <               Phunc * dlm_i(l,m) * dctidux
637 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
638 <               Phunc * dlm_i(l,m) * dctiduy
639 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
640 <               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 648 | 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 689 | Line 706 | contains  
706               dPhuncdUy = coeff*(spi * dum_i(m-1)*dcpiduy + dspiduy *um_i(m-1))
707               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
708            endif
709 <          !write(*,*) eps_i, plm_i(l,m), Phunc
710 <          eps_i = eps_i + plm_i(l,m)*Phunc
709 >
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 757 | Line 774 | contains  
774         xj2 = xj*xj
775         yj2 = yj*yj
776         zj2 = zj*zj
760 –      
761 –       projj = sqrt(xj2 + yj2)
762 –       projj3 = projj*projj*projj
763 –      
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
777 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
778 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
818 >       dcpjduz = 0.0d0
819        
820         spj = yj / projj
781 –       dspjdx = - xj * yj / projj3
782 –       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)
825 <      
826 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
827 <            ShapeMap%Shapes(st2)%bigM, LMAX, &
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 Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
# Line 840 | 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 882 | 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 924 | 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 977 | Line 1018 | contains  
1018      dsduxj = 0.5*dsjdux
1019      dsduyj = 0.5*dsjduy
1020      dsduzj = 0.5*dsjduz
1021 <    !write(*,*) eps_i, eps_j
1021 >
1022      eps = sqrt(eps_i * eps_j)
1023  
1024      depsdxi = eps_j * depsidx / (2.0d0 * eps)
# Line 994 | 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 1018 | 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 1044 | Line 1098 | contains  
1098      ! do the torques first since they are easy:
1099      ! remember that these are still in the body fixed axes
1100  
1047 –    txi = dvduxi * sw
1048 –    tyi = dvduyi * sw
1049 –    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 1115 | 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 1152 | 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)        
# Line 1311 | Line 1381 | contains  
1381         DY0 = DY1
1382         DY1 = DYN
1383      end DO
1384 +
1385 +
1386      RETURN
1387      
1388    end subroutine Orthogonal_Polynomial

Diff Legend

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