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 2211 by chrisfen, Thu Apr 21 14:12:19 2005 UTC vs.
Revision 2277 by chrisfen, Fri Aug 26 21:30:41 2005 UTC

# Line 69 | Line 69 | module shapes
69    public :: newShapeType
70    public :: complete_Shape_FF
71    public :: destroyShapeTypes
72 +  public :: getShapeCut
73  
74    type, private :: Shape
75       integer :: atid
# Line 97 | Line 98 | module shapes
98    type, private :: ShapeList
99       integer :: n_shapes = 0
100       integer :: currentShape = 0
101 <     type (Shape), pointer :: Shapes(:)      => null()
102 <     integer, pointer      :: atidToShape(:) => null()
101 >     type(Shape), pointer :: Shapes(:)      => null()
102 >     integer, pointer     :: atidToShape(:) => null()
103    end type ShapeList
104 <
104 >  
105    type(ShapeList), save :: ShapeMap
106 <
106 >  
107    integer :: lmax
108 <
108 >  
109   contains  
110 <
110 >  
111    subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
112         ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
113         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
114         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
115         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
116         c_ident, status)
117 <
117 >    
118      integer :: nContactFuncs
119      integer :: nRangeFuncs
120      integer :: nStrengthFuncs
# Line 154 | Line 155 | contains  
155         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
156  
157         ntypes = getSize(atypes)
158 <
158 >
159         allocate(ShapeMap%atidToShape(ntypes))
160      end if
161  
# Line 167 | Line 168 | contains  
168         status = -1
169         return
170      endif
170 –    
171 –    myATID = getFirstMatchingElement(atypes, 'c_ident', c_ident)
172 – !    write(*,*) 'myATID= ', myATID, ' c_ident = ', c_ident
171  
172 +    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
173 +
174      ShapeMap%atidToShape(myATID)                     = current
175      ShapeMap%Shapes(current)%atid                    = myATID
176      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
# Line 360 | Line 360 | contains  
360      
361         myATID = getFirstMatchingElement(atypes, 'c_ident', i)
362         call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
363 < !       write(*,*) 'is_LJ = ', thisProperty, ' for atid = ', myATID
364 <      
363 >        
364         if (thisProperty) then
366 –                
365            ShapeMap%currentShape = ShapeMap%currentShape + 1
366            current = ShapeMap%currentShape
367  
# Line 387 | Line 385 | contains  
385  
386    end subroutine complete_Shape_FF
387  
388 +  function getShapeCut(atomID) result(cutValue)
389 +    integer, intent(in) :: atomID
390 +    real(kind=dp) :: cutValue, whoopdedoo
391 +
392 +    !! this is just a placeholder for a cutoff value, hopefully we'll
393 +    !! develop a method to calculate a sensible value
394 +    whoopdedoo = 9.0_dp
395 +
396 +    cutValue = whoopdedoo
397 +
398 +  end function getShapeCut
399 +
400    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
401         pot, A, f, t, do_pot)
402  
# Line 487 | Line 497 | contains  
497      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
498      real (kind=dp) :: fxradial, fyradial, fzradial
499  
500 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
501 +
502      real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
503      real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
504      real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
# Line 505 | Line 517 | contains  
517      drdxi = -d(1) / rij
518      drdyi = -d(2) / rij
519      drdzi = -d(3) / rij
520 +    drduxi = 0.0d0
521 +    drduyi = 0.0d0
522 +    drduzi = 0.0d0
523  
524      drdxj = d(1) / rij
525      drdyj = d(2) / rij
526      drdzj = d(3) / rij
527 +    drduxj = 0.0d0
528 +    drduyj = 0.0d0
529 +    drduzj = 0.0d0
530  
531      ! find the atom type id (atid) for each atom:
532   #ifdef IS_MPI
# Line 567 | Line 585 | contains  
585         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
586  
587   #endif
588 <
588 >       xihat = xi / rij
589 >       yihat = yi / rij
590 >       zihat = zi / rij
591         xi2 = xi*xi
592         yi2 = yi*yi
593         zi2 = zi*zi            
# Line 579 | Line 599 | contains  
599         dctidx = - zi * xi / r3
600         dctidy = - zi * yi / r3
601         dctidz = 1.0d0 / rij - zi2 / r3
602 <       dctidux = - (zi * xi2) / r3
603 <       dctiduy = - (zi * yi2) / r3
604 <       dctiduz = zi / rij - (zi2 * zi) / r3
602 >       dctidux = yi / rij ! - (zi * xi2) / r3
603 >       dctiduy = -xi / rij !- (zi * yi2) / r3
604 >       dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3
605  
606         ! this is an attempt to try to truncate the singularity when
607         ! sin(theta) is near 0.0:
# Line 660 | Line 680 | contains  
680               dPhuncdX = coeff * dtm_i(m) * dcpidx
681               dPhuncdY = coeff * dtm_i(m) * dcpidy
682               dPhuncdZ = coeff * dtm_i(m) * dcpidz
683 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
683 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
684               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
685               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
686            else
# Line 674 | Line 694 | contains  
694            endif
695  
696            sigma_i = sigma_i + plm_i(m,l)*Phunc
697 <          write(*,*) 'dsigmaidux = ', dsigmaidux
698 <          write(*,*) 'Phunc = ', Phunc
697 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux
698 > !!$          write(*,*) 'Phunc = ', Phunc
699            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
700                 Phunc * dlm_i(m,l) * dctidx
701            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
# Line 688 | Line 708 | contains  
708                 Phunc * dlm_i(m,l) * dctiduy
709            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
710                 Phunc * dlm_i(m,l) * dctiduz
711 <          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
712 <                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
713 <                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
711 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
712 > !!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
713 > !!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
714         end do
715  
716         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 704 | Line 724 | contains  
724               dPhuncdX = coeff * dtm_i(m) * dcpidx
725               dPhuncdY = coeff * dtm_i(m) * dcpidy
726               dPhuncdZ = coeff * dtm_i(m) * dcpidz
727 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
727 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
728               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
729               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
730            else
# Line 746 | Line 766 | contains  
766               dPhuncdX = coeff * dtm_i(m) * dcpidx
767               dPhuncdY = coeff * dtm_i(m) * dcpidy
768               dPhuncdZ = coeff * dtm_i(m) * dcpidz
769 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
769 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
770               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
771               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
772            else
# Line 823 | Line 843 | contains  
843         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
844   #endif
845  
846 +       xjhat = xj / rij
847 +       yjhat = yj / rij
848 +       zjhat = zj / rij
849         xj2 = xj*xj
850         yj2 = yj*yj
851         zj2 = zj*zj
# Line 834 | Line 857 | contains  
857         dctjdx = - zj * xj / r3
858         dctjdy = - zj * yj / r3
859         dctjdz = 1.0d0 / rij - zj2 / r3
860 <       dctjdux = - (zi * xj2) / r3
861 <       dctjduy = - (zj * yj2) / r3
862 <       dctjduz = zj / rij - (zj2 * zj) / r3
860 >       dctjdux = yj / rij !- (zi * xj2) / r3
861 >       dctjduy = -xj / rij !- (zj * yj2) / r3
862 >       dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3
863  
864         ! this is an attempt to try to truncate the singularity when
865         ! sin(theta) is near 0.0:
# Line 918 | Line 941 | contains  
941               dPhuncdX = coeff * dtm_j(m) * dcpjdx
942               dPhuncdY = coeff * dtm_j(m) * dcpjdy
943               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
944 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
944 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
945               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
946               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
947            else
# Line 960 | Line 983 | contains  
983               dPhuncdX = coeff * dtm_j(m) * dcpjdx
984               dPhuncdY = coeff * dtm_j(m) * dcpjdy
985               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
986 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
986 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
987               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
988               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
989            else
# Line 1072 | Line 1095 | contains  
1095      dsduzj = 0.5*dsjduz
1096  
1097      eps = sqrt(eps_i * eps_j)
1098 <    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1099 <    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1100 < !    write(*,*) sigma_i, ' is sigma i; ', s_i, ' is s i; ', eps_i, ' is eps i'
1098 > !!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1099 > !!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1100 > !!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1101      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1102      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1103      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1089 | Line 1112 | contains  
1112      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1113      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1114  
1115 < !    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 < !    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1115 > !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1116  
1117 < !    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1118 < !    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1117 > !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1118 > !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1119   !!$
1120   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1121  
1122      rtdenom = rij-sigma+s
1123      rt = s / rtdenom
1124  
1125 <    drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1126 <    drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1127 <    drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1128 <    drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1129 <    write(*,*) dsduxi, ' is dsduxi; ', drduxi, ' is drduxi; ', dsigmaduxi, &
1130 <               ' is dsigmaduxi; ', dsduxi, ' is dsduxi'
1131 <    drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1132 <    drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1133 <    drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1134 <    drtdyj = (dsdyj + rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1135 <    drtdzj = (dsdzj + rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1136 <    drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1115 <    drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1116 <    drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1125 >    drtdxi = (dsdxi - rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1126 >    drtdyi = (dsdyi - rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1127 >    drtdzi = (dsdzi - rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1128 >    drtduxi = (dsduxi - rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1129 >    drtduyi = (dsduyi - rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1130 >    drtduzi = (dsduzi - rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1131 >    drtdxj = (dsdxj - rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1132 >    drtdyj = (dsdyj - rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1133 >    drtdzj = (dsdzj - rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1134 >    drtduxj = (dsduxj - rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1135 >    drtduyj = (dsduyj - rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1136 >    drtduzj = (dsduzj - rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1137  
1138 < !    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1139 < !    write(*,*) 'drtdu_i = ', drtduxi, drtduyi, drtduzi
1138 > !!$    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1139 > !!$    write(*,*) 'drtdu_j = ', drtduxj, drtduyj, drtduzj
1140  
1141      rt2 = rt*rt
1142      rt3 = rt2*rt
# Line 1153 | Line 1173 | contains  
1173      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1174      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1175      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1176 <    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1176 > !!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1177      ! do the torques first since they are easy:
1178      ! remember that these are still in the body fixed axes
1179  
1180 +    txi = 0.0d0
1181 +    tyi = 0.0d0
1182 +    tzi = 0.0d0
1183  
1184 < !!$    write(*,*) 'sw = ', sw
1185 < !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1186 < !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1164 < !!$
1165 < !    txi =  (dvduzi - dvduyi) * sw
1166 < !    tyi =  (dvduxi - dvduzi) * sw
1167 < !    tzi =  (dvduyi - dvduxi) * sw
1168 <
1169 < !    txj = (dvduzj - dvduyj) * sw
1170 < !    tyj = (dvduxj - dvduzj) * sw
1171 < !    tzj = (dvduyj - dvduxj) * sw
1184 >    txj = 0.0d0
1185 >    tyj = 0.0d0
1186 >    tzj = 0.0d0
1187  
1188 <    txi = dvduxi * sw
1189 <    tyi = dvduyi * sw
1190 <    tzi = dvduzi * sw
1188 >    txi = (dvduyi - dvduzi) * sw
1189 >    tyi = (dvduzi - dvduxi) * sw
1190 >    tzi = (dvduxi - dvduyi) * sw
1191  
1192 <    txj = dvduxj * sw
1193 <    tyj = dvduyj * sw
1194 <    tzj = dvduzj * sw
1192 >    txj = (dvduyj - dvduzj) * sw
1193 >    tyj = (dvduzj - dvduxj) * sw
1194 >    tzj = (dvduxj - dvduyj) * sw
1195  
1196 + !!$    txi = dvduxi * sw
1197 + !!$    tyi = dvduyi * sw
1198 + !!$    tzi = dvduzi * sw
1199 + !!$
1200 + !!$    txj = dvduxj * sw
1201 + !!$    tyj = dvduyj * sw
1202 + !!$    tzj = dvduzj * sw
1203 +
1204      write(*,*) 't1 = ', txi, tyi, tzi
1205      write(*,*) 't2 = ', txj, tyj, tzj
1206  
# Line 1205 | Line 1228 | contains  
1228      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1229      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1230      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1231 +
1232   #endif
1233      ! Now, on to the forces:
1234  
1235      ! first rotate the i terms back into the lab frame:
1236  
1237 <    fxi = dvdxi * sw
1238 <    fyi = dvdyi * sw
1239 <    fzi = dvdzi * sw
1237 >    fxi = -dvdxi * sw
1238 >    fyi = -dvdyi * sw
1239 >    fzi = -dvdzi * sw
1240  
1241 <    fxj = dvdxj * sw
1242 <    fyj = dvdyj * sw
1243 <    fzj = dvdzj * sw
1241 >    fxj = -dvdxj * sw
1242 >    fyj = -dvdyj * sw
1243 >    fzj = -dvdzj * sw
1244  
1245  
1246   #ifdef IS_MPI
# Line 1245 | Line 1269 | contains  
1269      fyji = -fyjj
1270      fzji = -fzjj
1271  
1272 <    fxradial = 0.5_dp * (fxii + fxji)
1273 <    fyradial = 0.5_dp * (fyii + fyji)
1274 <    fzradial = 0.5_dp * (fzii + fzji)
1275 <    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1272 >    fxradial = (fxii + fxji)
1273 >    fyradial = (fyii + fyji)
1274 >    fzradial = (fzii + fzji)
1275 > !!$    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1276   #ifdef IS_MPI
1277      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1278      f_Row(2,atom1) = f_Row(2,atom1) + fyradial

Diff Legend

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