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

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1707 by gezelter, Thu Nov 4 16:20:37 2004 UTC vs.
Revision 1717 by chrisfen, Fri Nov 5 21:04:33 2004 UTC

# Line 348 | 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 359 | 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 446 | Line 447 | contains  
447         call handleError("calc_shape", "NO SHAPEMAP!!!!")
448         return      
449      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)
450      
451      !! We assume that the rotation matrices have already been calculated
452      !! and placed in the A array.
# Line 535 | Line 529 | contains  
529         dctidx = - zi * xi / r3
530         dctidy = - zi * yi / r3
531         dctidz = 1.0d0 / rij - zi2 / r3
532 <       dctidux =  yi / rij + (zi * yi) / r3
533 <       dctiduy = -xi / rij - (zi * xi) / r3
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:
# Line 547 | Line 541 | contains  
541            proji = sqrt(rij * 1.0d-12)
542            dcpidx = 1.0d0 / proji
543            dcpidy = 0.0d0
544 <          dcpidux = 0.0d0
545 <          dcpiduy = zi / proji
544 >          dcpidux = xi / proji
545 >          dcpiduy = 0.0d0
546            dspidx = 0.0d0
547            dspidy = 1.0d0 / proji
548 <          dspidux = -zi / proji
549 <          dspiduy = 0.0d0
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 * yi * zi / proji3
556 <          dcpiduy = zi / proji - xi2 * zi / 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 = -zi / proji + yi2 * zi / proji3
560 <          dspiduy = - xi * yi * zi / proji3
559 >          dspidux = - (yi * xi2) / proji3
560 >          dspiduy = yi / proji - (yi2 * yi) / proji3
561         endif
562        
563         cpi = xi / proji
564         dcpidz = 0.0d0
565 <       dcpiduz = -yi / proji
565 >       dcpiduz = 0.0d0
566        
567         spi = yi / proji
568         dspidz = 0.0d0
569 <       dspiduz = xi / proji
576 <       write(*,*) 'before lmloop', cpi, dcpidx, dcpidux
569 >       dspiduz = 0.0d0
570  
571         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
572              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 654 | Line 647 | contains  
647            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
648            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
649            
657 –         write(*,*) 'in lm loop a', coeff, dtm_i(m), dcpidx
658 –
659 –
650            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
651               Phunc = coeff * tm_i(m)
652               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 676 | Line 666 | contains  
666            endif
667  
668            s_i = s_i + plm_i(m,l)*Phunc
669 <          
680 <
681 <          write(*,*) 'in lm loop ', dsidx, plm_i(m,l), dPhuncdX, Phunc, dlm_i(m,l), dctidx
669 >          
670            dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
671                 Phunc * dlm_i(m,l) * dctidx
672            dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
# Line 794 | Line 782 | contains  
782         dctjdx = - zj * xj / r3
783         dctjdy = - zj * yj / r3
784         dctjdz = 1.0d0 / rij - zj2 / r3
785 <       dctjdux =  yj / rij + (zj * yj) / r3
786 <       dctjduy = -xj / rij - (zj * xj) / r3
787 <       dctjduz = 0.0d0
785 >       dctjdux = - (zi * xj2) / r3
786 >       dctjduy = - (zj * yj2) / r3
787 >       dctjduz = zj / rij - (zj2 * zj) / r3
788        
789         ! this is an attempt to try to truncate the singularity when
790         ! sin(theta) is near 0.0:
# Line 806 | Line 794 | contains  
794            projj = sqrt(rij * 1.0d-12)
795            dcpjdx = 1.0d0 / projj
796            dcpjdy = 0.0d0
797 <          dcpjdux = 0.0d0
798 <          dcpjduy = zj / projj
797 >          dcpjdux = xj / projj
798 >          dcpjduy = 0.0d0
799            dspjdx = 0.0d0
800            dspjdy = 1.0d0 / projj
801 <          dspjdux = -zj / projj
802 <          dspjduy = 0.0d0
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 * yj * zj / projj3
809 <          dcpjduy = zj / projj - xj2 * zj / 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 = -zj / projj + yj2 * zj / projj3
813 <          dspjduy = - xj * yj * zj / 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 <       dcpjduz = -yj / projj
818 >       dcpjduz = 0.0d0
819        
820         spj = yj / projj
821         dspjdz = 0.0d0
822 <       dspjduz = xj / projj
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)
# Line 972 | Line 963 | contains  
963               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
964            endif
965  
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(m,l)*dPhuncdX + &
# Line 1028 | Line 1021 | contains  
1021  
1022      eps = sqrt(eps_i * eps_j)
1023  
1031 –    write(*,*) 'sigma, s, eps = ', sigma, s, eps
1032 –
1024      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1025      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1026      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1044 | Line 1035 | contains  
1035      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1036      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1037      
1038 <    rtdenom = rij-sigma+s
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 <    write(*,*) 'rtdenom = ', rtdenom, ' sw = ', sw
1046 >    rtdenom = rij-sigma+s
1047      rt = s / rtdenom
1048  
1052 –    write(*,*) 'john' , dsdxi, rt, drdxi, dsigmadxi, rtdenom
1053 –    write(*,*) 'bigboot', dsduzj, rt, drduzj, dsigmaduzj, rtdenom
1054 –
1055 –
1049      drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1050      drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1051      drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
# Line 1074 | 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
1085 –    
1086 –    write(*,*) 'drtdxi = ', drtdxi, drtdyi
1087 –    write(*,*) 'depsdxi = ', depsdxi, depsdyi
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
1086      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1093 | Line 1088 | contains  
1088      dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1089      dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1090  
1096 –    write(*,*) 'drtduzj = ', drtduzj, depsduzj
1097 –
1091      dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1092      dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1093      dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
# Line 1105 | Line 1098 | contains  
1098      ! do the torques first since they are easy:
1099      ! remember that these are still in the body fixed axes
1100  
1108 –    write(*,*) 'dvdx = ', dvdxi, dvdyi, dvdzi
1109 –    write(*,*) 'dvdx = ', dvdxj, dvdyj, dvdzj
1110 –    write(*,*) 'dvdu = ', dvduxi, dvduyi, dvduzi
1111 –    write(*,*) 'dvdu = ', dvduxj, dvduyj, dvduzj
1101  
1102 <    txi = dvduxi * sw
1103 <    tyi = dvduyi * sw
1104 <    tzi = dvduzi * 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 = dvduxj * sw
1111 <    tyj = dvduyj * sw
1112 <    tzj = dvduzj * sw
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 1181 | 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 1218 | Line 1222 | contains  
1222         fpair(3) = fpair(3) + fzradial
1223        
1224      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)
1225  
1227 –
1226    end subroutine do_shape_pair
1227      
1228    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        

Diff Legend

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