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

Comparing trunk/OOPSE-2.0/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 2213 by chrisfen, Thu Apr 21 14:12:19 2005 UTC vs.
Revision 2214 by chrisfen, Wed Apr 27 20:14:03 2005 UTC

# Line 97 | Line 97 | module shapes
97    type, private :: ShapeList
98       integer :: n_shapes = 0
99       integer :: currentShape = 0
100 <     type (Shape), pointer :: Shapes(:)      => null()
101 <     integer, pointer      :: atidToShape(:) => null()
100 >     type(Shape), pointer :: Shapes(:)      => null()
101 >     integer, pointer     :: atidToShape(:) => null()
102    end type ShapeList
103 <
103 >  
104    type(ShapeList), save :: ShapeMap
105 <
105 >  
106    integer :: lmax
107 <
107 >  
108   contains  
109 <
109 >  
110    subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
111         ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115         c_ident, status)
116 <
116 >    
117      integer :: nContactFuncs
118      integer :: nRangeFuncs
119      integer :: nStrengthFuncs
# Line 154 | Line 154 | contains  
154         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
155  
156         ntypes = getSize(atypes)
157 <
157 >
158         allocate(ShapeMap%atidToShape(ntypes))
159      end if
160  
# Line 167 | Line 167 | contains  
167         status = -1
168         return
169      endif
170 –    
171 –    myATID = getFirstMatchingElement(atypes, 'c_ident', c_ident)
172 – !    write(*,*) 'myATID= ', myATID, ' c_ident = ', c_ident
170  
171 +    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
172 +
173      ShapeMap%atidToShape(myATID)                     = current
174      ShapeMap%Shapes(current)%atid                    = myATID
175      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
# Line 360 | Line 359 | contains  
359      
360         myATID = getFirstMatchingElement(atypes, 'c_ident', i)
361         call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
362 < !       write(*,*) 'is_LJ = ', thisProperty, ' for atid = ', myATID
364 <      
362 >        
363         if (thisProperty) then
366 –                
364            ShapeMap%currentShape = ShapeMap%currentShape + 1
365            current = ShapeMap%currentShape
366  
# Line 487 | Line 484 | contains  
484      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
485      real (kind=dp) :: fxradial, fyradial, fzradial
486  
487 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
488 +
489      real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
490      real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
491      real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
# Line 505 | Line 504 | contains  
504      drdxi = -d(1) / rij
505      drdyi = -d(2) / rij
506      drdzi = -d(3) / rij
507 +    drduxi = 0.0d0
508 +    drduyi = 0.0d0
509 +    drduzi = 0.0d0
510  
511      drdxj = d(1) / rij
512      drdyj = d(2) / rij
513      drdzj = d(3) / rij
514 +    drduxj = 0.0d0
515 +    drduyj = 0.0d0
516 +    drduzj = 0.0d0
517  
518      ! find the atom type id (atid) for each atom:
519   #ifdef IS_MPI
# Line 567 | Line 572 | contains  
572         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
573  
574   #endif
575 <
575 >       xihat = xi / rij
576 >       yihat = yi / rij
577 >       zihat = zi / rij
578         xi2 = xi*xi
579         yi2 = yi*yi
580         zi2 = zi*zi            
# Line 579 | Line 586 | contains  
586         dctidx = - zi * xi / r3
587         dctidy = - zi * yi / r3
588         dctidz = 1.0d0 / rij - zi2 / r3
589 <       dctidux = - (zi * xi2) / r3
590 <       dctiduy = - (zi * yi2) / r3
591 <       dctiduz = zi / rij - (zi2 * zi) / r3
589 >       dctidux = yi / rij ! - (zi * xi2) / r3
590 >       dctiduy = -xi / rij !- (zi * yi2) / r3
591 >       dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3
592  
593         ! this is an attempt to try to truncate the singularity when
594         ! sin(theta) is near 0.0:
# Line 660 | Line 667 | contains  
667               dPhuncdX = coeff * dtm_i(m) * dcpidx
668               dPhuncdY = coeff * dtm_i(m) * dcpidy
669               dPhuncdZ = coeff * dtm_i(m) * dcpidz
670 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
670 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
671               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
672               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
673            else
# Line 674 | Line 681 | contains  
681            endif
682  
683            sigma_i = sigma_i + plm_i(m,l)*Phunc
684 <          write(*,*) 'dsigmaidux = ', dsigmaidux
685 <          write(*,*) 'Phunc = ', Phunc
684 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux
685 > !!$          write(*,*) 'Phunc = ', Phunc
686            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
687                 Phunc * dlm_i(m,l) * dctidx
688            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
# Line 688 | Line 695 | contains  
695                 Phunc * dlm_i(m,l) * dctiduy
696            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
697                 Phunc * dlm_i(m,l) * dctiduz
698 <          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
699 <                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
700 <                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
698 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
699 > !!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
700 > !!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
701         end do
702  
703         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 704 | Line 711 | contains  
711               dPhuncdX = coeff * dtm_i(m) * dcpidx
712               dPhuncdY = coeff * dtm_i(m) * dcpidy
713               dPhuncdZ = coeff * dtm_i(m) * dcpidz
714 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
714 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
715               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
716               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
717            else
# Line 746 | Line 753 | contains  
753               dPhuncdX = coeff * dtm_i(m) * dcpidx
754               dPhuncdY = coeff * dtm_i(m) * dcpidy
755               dPhuncdZ = coeff * dtm_i(m) * dcpidz
756 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
756 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
757               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
758               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
759            else
# Line 823 | Line 830 | contains  
830         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
831   #endif
832  
833 +       xjhat = xj / rij
834 +       yjhat = yj / rij
835 +       zjhat = zj / rij
836         xj2 = xj*xj
837         yj2 = yj*yj
838         zj2 = zj*zj
# Line 834 | Line 844 | contains  
844         dctjdx = - zj * xj / r3
845         dctjdy = - zj * yj / r3
846         dctjdz = 1.0d0 / rij - zj2 / r3
847 <       dctjdux = - (zi * xj2) / r3
848 <       dctjduy = - (zj * yj2) / r3
849 <       dctjduz = zj / rij - (zj2 * zj) / r3
847 >       dctjdux = yj / rij !- (zi * xj2) / r3
848 >       dctjduy = -xj / rij !- (zj * yj2) / r3
849 >       dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3
850  
851         ! this is an attempt to try to truncate the singularity when
852         ! sin(theta) is near 0.0:
# Line 918 | Line 928 | contains  
928               dPhuncdX = coeff * dtm_j(m) * dcpjdx
929               dPhuncdY = coeff * dtm_j(m) * dcpjdy
930               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
931 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
931 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
932               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
933               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
934            else
# Line 960 | Line 970 | contains  
970               dPhuncdX = coeff * dtm_j(m) * dcpjdx
971               dPhuncdY = coeff * dtm_j(m) * dcpjdy
972               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
973 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
973 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
974               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
975               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
976            else
# Line 1072 | Line 1082 | contains  
1082      dsduzj = 0.5*dsjduz
1083  
1084      eps = sqrt(eps_i * eps_j)
1085 <    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1086 <    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1087 < !    write(*,*) sigma_i, ' is sigma i; ', s_i, ' is s i; ', eps_i, ' is eps i'
1085 > !!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1086 > !!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1087 > !!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1088      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1089      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1090      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1089 | Line 1099 | contains  
1099      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1100      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1101  
1102 < !    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 < !    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1102 > !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1103  
1104 < !    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1105 < !    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1104 > !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1105 > !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1106   !!$
1107   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1108  
1109      rtdenom = rij-sigma+s
1110      rt = s / rtdenom
1111  
1112 <    drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1113 <    drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1114 <    drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1115 <    drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1116 <    write(*,*) dsduxi, ' is dsduxi; ', drduxi, ' is drduxi; ', dsigmaduxi, &
1117 <               ' is dsigmaduxi; ', dsduxi, ' is dsduxi'
1118 <    drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1119 <    drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1120 <    drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1121 <    drtdyj = (dsdyj + rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1122 <    drtdzj = (dsdzj + rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1123 <    drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1115 <    drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1116 <    drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1112 >    drtdxi = (dsdxi - rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1113 >    drtdyi = (dsdyi - rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1114 >    drtdzi = (dsdzi - rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1115 >    drtduxi = (dsduxi - rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1116 >    drtduyi = (dsduyi - rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1117 >    drtduzi = (dsduzi - rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1118 >    drtdxj = (dsdxj - rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1119 >    drtdyj = (dsdyj - rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1120 >    drtdzj = (dsdzj - rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1121 >    drtduxj = (dsduxj - rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1122 >    drtduyj = (dsduyj - rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1123 >    drtduzj = (dsduzj - rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1124  
1125 < !    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1126 < !    write(*,*) 'drtdu_i = ', drtduxi, drtduyi, drtduzi
1125 > !!$    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1126 > !!$    write(*,*) 'drtdu_j = ', drtduxj, drtduyj, drtduzj
1127  
1128      rt2 = rt*rt
1129      rt3 = rt2*rt
# Line 1153 | Line 1160 | contains  
1160      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1161      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1162      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1163 <    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1163 > !!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1164      ! do the torques first since they are easy:
1165      ! remember that these are still in the body fixed axes
1166  
1167 +    txi = 0.0d0
1168 +    tyi = 0.0d0
1169 +    tzi = 0.0d0
1170  
1171 < !!$    write(*,*) 'sw = ', sw
1172 < !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1173 < !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1164 < !!$
1165 < !    txi =  (dvduzi - dvduyi) * sw
1166 < !    tyi =  (dvduxi - dvduzi) * sw
1167 < !    tzi =  (dvduyi - dvduxi) * sw
1171 >    txj = 0.0d0
1172 >    tyj = 0.0d0
1173 >    tzj = 0.0d0
1174  
1175 < !    txj = (dvduzj - dvduyj) * sw
1176 < !    tyj = (dvduxj - dvduzj) * sw
1177 < !    tzj = (dvduyj - dvduxj) * sw
1175 >    txi = (dvduyi - dvduzi) * sw
1176 >    tyi = (dvduzi - dvduxi) * sw
1177 >    tzi = (dvduxi - dvduyi) * sw
1178  
1179 <    txi = dvduxi * sw
1180 <    tyi = dvduyi * sw
1181 <    tzi = dvduzi * sw
1179 >    txj = (dvduyj - dvduzj) * sw
1180 >    tyj = (dvduzj - dvduxj) * sw
1181 >    tzj = (dvduxj - dvduyj) * sw
1182  
1183 <    txj = dvduxj * sw
1184 <    tyj = dvduyj * sw
1185 <    tzj = dvduzj * sw
1183 > !!$    txi = dvduxi * sw
1184 > !!$    tyi = dvduyi * sw
1185 > !!$    tzi = dvduzi * sw
1186 > !!$
1187 > !!$    txj = dvduxj * sw
1188 > !!$    tyj = dvduyj * sw
1189 > !!$    tzj = dvduzj * sw
1190  
1191      write(*,*) 't1 = ', txi, tyi, tzi
1192      write(*,*) 't2 = ', txj, tyj, tzj
# Line 1205 | Line 1215 | contains  
1215      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1216      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1217      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1218 +
1219   #endif
1220      ! Now, on to the forces:
1221  
1222      ! first rotate the i terms back into the lab frame:
1223  
1224 <    fxi = dvdxi * sw
1225 <    fyi = dvdyi * sw
1226 <    fzi = dvdzi * sw
1224 >    fxi = -dvdxi * sw
1225 >    fyi = -dvdyi * sw
1226 >    fzi = -dvdzi * sw
1227  
1228 <    fxj = dvdxj * sw
1229 <    fyj = dvdyj * sw
1230 <    fzj = dvdzj * sw
1228 >    fxj = -dvdxj * sw
1229 >    fyj = -dvdyj * sw
1230 >    fzj = -dvdzj * sw
1231  
1232  
1233   #ifdef IS_MPI
# Line 1245 | Line 1256 | contains  
1256      fyji = -fyjj
1257      fzji = -fzjj
1258  
1259 <    fxradial = 0.5_dp * (fxii + fxji)
1260 <    fyradial = 0.5_dp * (fyii + fyji)
1261 <    fzradial = 0.5_dp * (fzii + fzji)
1262 <    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1259 >    fxradial = (fxii + fxji)
1260 >    fyradial = (fyii + fyji)
1261 >    fzradial = (fzii + fzji)
1262 > !!$    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1263   #ifdef IS_MPI
1264      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1265      f_Row(2,atom1) = f_Row(2,atom1) + fyradial

Diff Legend

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