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 1707 by gezelter, Thu Nov 4 16:20:37 2004 UTC vs.
Revision 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC

# Line 1 | Line 1
1 + !!
2 + !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + !!
4 + !! The University of Notre Dame grants you ("Licensee") a
5 + !! non-exclusive, royalty free, license to use, modify and
6 + !! redistribute this software in source and binary code form, provided
7 + !! that the following conditions are met:
8 + !!
9 + !! 1. Acknowledgement of the program authors must be made in any
10 + !!    publication of scientific results based in part on use of the
11 + !!    program.  An acceptable form of acknowledgement is citation of
12 + !!    the article in which the program was described (Matthew
13 + !!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + !!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + !!    Parallel Simulation Engine for Molecular Dynamics,"
16 + !!    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + !!
18 + !! 2. Redistributions of source code must retain the above copyright
19 + !!    notice, this list of conditions and the following disclaimer.
20 + !!
21 + !! 3. Redistributions in binary form must reproduce the above copyright
22 + !!    notice, this list of conditions and the following disclaimer in the
23 + !!    documentation and/or other materials provided with the
24 + !!    distribution.
25 + !!
26 + !! This software is provided "AS IS," without a warranty of any
27 + !! kind. All express or implied conditions, representations and
28 + !! warranties, including any implied warranty of merchantability,
29 + !! fitness for a particular purpose or non-infringement, are hereby
30 + !! excluded.  The University of Notre Dame and its licensors shall not
31 + !! be liable for any damages suffered by licensee as a result of
32 + !! using, modifying or distributing the software or its
33 + !! derivatives. In no event will the University of Notre Dame or its
34 + !! licensors be liable for any lost revenue, profit or data, or for
35 + !! direct, indirect, special, consequential, incidental or punitive
36 + !! damages, however caused and regardless of the theory of liability,
37 + !! arising out of the use of or inability to use software, even if the
38 + !! University of Notre Dame has been advised of the possibility of
39 + !! such damages.
40 + !!
41 +
42 +
43   module shapes
44  
45    use force_globals
# Line 348 | Line 390 | contains  
390      real (kind=dp), intent(inout) :: rij, r2
391      real (kind=dp), dimension(3), intent(in) :: d
392      real (kind=dp), dimension(3), intent(inout) :: fpair
393 <    real (kind=dp) :: pot, vpair, sw
393 >    real (kind=dp) :: pot, vpair, sw, dswdr
394      real (kind=dp), dimension(9,nLocal) :: A
395      real (kind=dp), dimension(3,nLocal) :: f
396      real (kind=dp), dimension(3,nLocal) :: t
# Line 359 | Line 401 | contains  
401      integer :: l, m, lm, id1, id2, localError, function_type
402      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
403      real (kind=dp) :: coeff
404 +    real (kind=dp) :: pot_temp
405  
406      real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
407      real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
# Line 446 | Line 489 | contains  
489         call handleError("calc_shape", "NO SHAPEMAP!!!!")
490         return      
491      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)
492      
493      !! We assume that the rotation matrices have already been calculated
494      !! and placed in the A array.
# Line 535 | Line 571 | contains  
571         dctidx = - zi * xi / r3
572         dctidy = - zi * yi / r3
573         dctidz = 1.0d0 / rij - zi2 / r3
574 <       dctidux =  yi / rij + (zi * yi) / r3
575 <       dctiduy = -xi / rij - (zi * xi) / r3
576 <       dctiduz = 0.0d0
574 >       dctidux = - (zi * xi2) / r3
575 >       dctiduy = - (zi * yi2) / r3
576 >       dctiduz = zi / rij - (zi2 * zi) / r3
577  
578         ! this is an attempt to try to truncate the singularity when
579         ! sin(theta) is near 0.0:
# Line 547 | Line 583 | contains  
583            proji = sqrt(rij * 1.0d-12)
584            dcpidx = 1.0d0 / proji
585            dcpidy = 0.0d0
586 <          dcpidux = 0.0d0
587 <          dcpiduy = zi / proji
586 >          dcpidux = xi / proji
587 >          dcpiduy = 0.0d0
588            dspidx = 0.0d0
589            dspidy = 1.0d0 / proji
590 <          dspidux = -zi / proji
591 <          dspiduy = 0.0d0
590 >          dspidux = 0.0d0
591 >          dspiduy = yi / proji
592         else
593            proji = sqrt(xi2 + yi2)
594            proji3 = proji*proji*proji
595            dcpidx = 1.0d0 / proji - xi2 / proji3
596            dcpidy = - xi * yi / proji3
597 <          dcpidux = xi * yi * zi / proji3
598 <          dcpiduy = zi / proji - xi2 * zi / proji3
597 >          dcpidux = xi / proji - (xi2 * xi) / proji3
598 >          dcpiduy = - (xi * yi2) / proji3
599            dspidx = - xi * yi / proji3
600            dspidy = 1.0d0 / proji - yi2 / proji3
601 <          dspidux = -zi / proji + yi2 * zi / proji3
602 <          dspiduy = - xi * yi * zi / proji3
601 >          dspidux = - (yi * xi2) / proji3
602 >          dspiduy = yi / proji - (yi2 * yi) / proji3
603         endif
604        
605         cpi = xi / proji
606         dcpidz = 0.0d0
607 <       dcpiduz = -yi / proji
607 >       dcpiduz = 0.0d0
608        
609         spi = yi / proji
610         dspidz = 0.0d0
611 <       dspiduz = xi / proji
576 <       write(*,*) 'before lmloop', cpi, dcpidx, dcpidux
611 >       dspiduz = 0.0d0
612  
613         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
614              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 654 | Line 689 | contains  
689            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
690            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
691            
657 –         write(*,*) 'in lm loop a', coeff, dtm_i(m), dcpidx
658 –
659 –
692            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
693               Phunc = coeff * tm_i(m)
694               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 676 | Line 708 | contains  
708            endif
709  
710            s_i = s_i + plm_i(m,l)*Phunc
711 <          
680 <
681 <          write(*,*) 'in lm loop ', dsidx, plm_i(m,l), dPhuncdX, Phunc, dlm_i(m,l), dctidx
711 >          
712            dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
713                 Phunc * dlm_i(m,l) * dctidx
714            dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
# Line 794 | Line 824 | contains  
824         dctjdx = - zj * xj / r3
825         dctjdy = - zj * yj / r3
826         dctjdz = 1.0d0 / rij - zj2 / r3
827 <       dctjdux =  yj / rij + (zj * yj) / r3
828 <       dctjduy = -xj / rij - (zj * xj) / r3
829 <       dctjduz = 0.0d0
827 >       dctjdux = - (zi * xj2) / r3
828 >       dctjduy = - (zj * yj2) / r3
829 >       dctjduz = zj / rij - (zj2 * zj) / r3
830        
831         ! this is an attempt to try to truncate the singularity when
832         ! sin(theta) is near 0.0:
# Line 806 | Line 836 | contains  
836            projj = sqrt(rij * 1.0d-12)
837            dcpjdx = 1.0d0 / projj
838            dcpjdy = 0.0d0
839 <          dcpjdux = 0.0d0
840 <          dcpjduy = zj / projj
839 >          dcpjdux = xj / projj
840 >          dcpjduy = 0.0d0
841            dspjdx = 0.0d0
842            dspjdy = 1.0d0 / projj
843 <          dspjdux = -zj / projj
844 <          dspjduy = 0.0d0
843 >          dspjdux = 0.0d0
844 >          dspjduy = yj / projj
845         else
846            projj = sqrt(xj2 + yj2)
847            projj3 = projj*projj*projj
848            dcpjdx = 1.0d0 / projj - xj2 / projj3
849            dcpjdy = - xj * yj / projj3
850 <          dcpjdux = xj * yj * zj / projj3
851 <          dcpjduy = zj / projj - xj2 * zj / projj3
850 >          dcpjdux = xj / projj - (xj2 * xj) / projj3
851 >          dcpjduy = - (xj * yj2) / projj3
852            dspjdx = - xj * yj / projj3
853            dspjdy = 1.0d0 / projj - yj2 / projj3
854 <          dspjdux = -zj / projj + yj2 * zj / projj3
855 <          dspjduy = - xj * yj * zj / projj3
854 >          dspjdux = - (yj * xj2) / projj3
855 >          dspjduy = yj / projj - (yj2 * yj) / projj3
856         endif
857  
858         cpj = xj / projj
859         dcpjdz = 0.0d0
860 <       dcpjduz = -yj / projj
860 >       dcpjduz = 0.0d0
861        
862         spj = yj / projj
863         dspjdz = 0.0d0
864 <       dspjduz = xj / projj
864 >       dspjduz = 0.0d0
865  
866 +
867 +       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
868 +       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
869         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
870              ShapeMap%Shapes(st2)%bigL, LMAX, &
871              plm_j, dlm_j)
# Line 972 | Line 1005 | contains  
1005               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1006            endif
1007  
1008 +          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1009 +
1010            eps_j = eps_j + plm_j(m,l)*Phunc
1011            
1012            depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
# Line 1028 | Line 1063 | contains  
1063  
1064      eps = sqrt(eps_i * eps_j)
1065  
1031 –    write(*,*) 'sigma, s, eps = ', sigma, s, eps
1032 –
1066      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1067      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1068      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1044 | Line 1077 | contains  
1077      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1078      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1079      
1080 <    rtdenom = rij-sigma+s
1080 > !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1081 > !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1082 > !!$
1083 > !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1084 > !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1085 > !!$
1086 > !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1087  
1088 <    write(*,*) 'rtdenom = ', rtdenom, ' sw = ', sw
1088 >    rtdenom = rij-sigma+s
1089      rt = s / rtdenom
1090  
1052 –    write(*,*) 'john' , dsdxi, rt, drdxi, dsigmadxi, rtdenom
1053 –    write(*,*) 'bigboot', dsduzj, rt, drduzj, dsigmaduzj, rtdenom
1054 –
1055 –
1091      drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1092      drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1093      drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
# Line 1074 | Line 1109 | contains  
1109      rt12 = rt6*rt6
1110      rt126 = rt12 - rt6
1111  
1112 +    pot_temp = 4.0d0 * eps * rt126
1113 +
1114 +    vpair = vpair + pot_temp
1115      if (do_pot) then
1116   #ifdef IS_MPI
1117 <       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1118 <       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1117 >       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1118 >       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1119   #else
1120 <       pot = pot + 4.0d0*eps*rt126*sw
1120 >       pot = pot + pot_temp*sw
1121   #endif
1122      endif
1085 –    
1086 –    write(*,*) 'drtdxi = ', drtdxi, drtdyi
1087 –    write(*,*) 'depsdxi = ', depsdxi, depsdyi
1123  
1124 + !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1125 +    
1126      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1127      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1128      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1093 | Line 1130 | contains  
1130      dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1131      dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1132  
1096 –    write(*,*) 'drtduzj = ', drtduzj, depsduzj
1097 –
1133      dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1134      dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1135      dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
# Line 1105 | Line 1140 | contains  
1140      ! do the torques first since they are easy:
1141      ! remember that these are still in the body fixed axes
1142  
1108 –    write(*,*) 'dvdx = ', dvdxi, dvdyi, dvdzi
1109 –    write(*,*) 'dvdx = ', dvdxj, dvdyj, dvdzj
1110 –    write(*,*) 'dvdu = ', dvduxi, dvduyi, dvduzi
1111 –    write(*,*) 'dvdu = ', dvduxj, dvduyj, dvduzj
1143  
1144 <    txi = dvduxi * sw
1145 <    tyi = dvduyi * sw
1146 <    tzi = dvduzi * sw
1144 > !!$    write(*,*) 'sw = ', sw
1145 > !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1146 > !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1147 > !!$
1148 >    txi =  (dvduzi - dvduyi) * sw
1149 >    tyi =  (dvduxi - dvduzi) * sw
1150 >    tzi =  (dvduyi - dvduxi) * sw
1151  
1152 <    txj = dvduxj * sw
1153 <    tyj = dvduyj * sw
1154 <    tzj = dvduzj * sw
1152 >    txj = (dvduzj - dvduyj) * sw
1153 >    tyj = (dvduxj - dvduzj) * sw
1154 >    tzj = (dvduyj - dvduxj) * sw
1155  
1156 + !!$    txi = -dvduxi * sw
1157 + !!$    tyi = -dvduyi * sw
1158 + !!$    tzi = -dvduzi * sw
1159 + !!$
1160 + !!$    txj = dvduxj * sw
1161 + !!$    tyj = dvduyj * sw
1162 + !!$    tzj = dvduzj * sw
1163 +
1164 +    write(*,*) 't1 = ', txi, tyi, tzi
1165 +    write(*,*) 't2 = ', txj, tyj, tzj
1166 +
1167      ! go back to lab frame using transpose of rotation matrix:
1168      
1169   #ifdef IS_MPI
# Line 1181 | Line 1227 | contains  
1227      fyji = -fyjj
1228      fzji = -fzjj
1229  
1230 <    fxradial = fxii + fxji
1231 <    fyradial = fyii + fyji
1232 <    fzradial = fzii + fzji
1230 >    fxradial = 0.5_dp * (fxii + fxji)
1231 >    fyradial = 0.5_dp * (fyii + fyji)
1232 >    fzradial = 0.5_dp * (fzii + fzji)
1233  
1234   #ifdef IS_MPI
1235      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
# Line 1218 | Line 1264 | contains  
1264         fpair(3) = fpair(3) + fzradial
1265        
1266      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)
1267  
1227 –
1268    end subroutine do_shape_pair
1269      
1270    SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
# Line 1390 | Line 1430 | end module shapes
1430    end subroutine Orthogonal_Polynomial
1431    
1432   end module shapes
1393 –
1394 – subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1395 –     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1396 –     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1397 –     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1398 –     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1399 –     myATID, status)
1400 –
1401 –  use definitions
1402 –  use shapes, only: newShapeType
1403 –  
1404 –  integer :: nContactFuncs
1405 –  integer :: nRangeFuncs
1406 –  integer :: nStrengthFuncs
1407 –  integer :: status
1408 –  integer :: myATID
1409 –  
1410 –  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1411 –  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1412 –  integer, dimension(nContactFuncs) :: ContactFunctionType        
1413 –  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1414 –  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1415 –  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1416 –  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1417 –  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1418 –  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1419 –  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1420 –  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1421 –  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1422 –  
1423 –  call newShapeType(nContactFuncs, ContactFuncLValue, &
1424 –       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1425 –       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1426 –       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1427 –       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1428 –       myATID, status)
1429 –
1430 –  return
1431 – end subroutine makeShape
1432 –
1433 – subroutine completeShapeFF(status)
1434 –
1435 –  use shapes, only: complete_Shape_FF
1436 –
1437 –  integer, intent(out)  :: status
1438 –  integer :: myStatus
1439 –
1440 –  myStatus = 0
1441 –
1442 –  call complete_Shape_FF(myStatus)
1443 –
1444 –  status = myStatus
1445 –
1446 –  return
1447 – end subroutine completeShapeFF
1448 –

Diff Legend

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