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 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2211 by chrisfen, Thu Apr 21 14:12:19 2005 UTC

# Line 112 | Line 112 | contains  
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115 <       myATID, status)
115 >       c_ident, status)
116  
117      integer :: nContactFuncs
118      integer :: nRangeFuncs
119      integer :: nStrengthFuncs
120      integer :: shape_ident
121      integer :: status
122 +    integer :: c_ident
123      integer :: myATID
124      integer :: bigL
125      integer :: bigM
# Line 154 | Line 155 | contains  
155  
156         ntypes = getSize(atypes)
157  
158 <       allocate(ShapeMap%atidToShape(0:ntypes))
158 >       allocate(ShapeMap%atidToShape(ntypes))
159      end if
160  
161      ShapeMap%currentShape = ShapeMap%currentShape + 1
# Line 166 | Line 167 | contains  
167         status = -1
168         return
169      endif
170 <
171 <    call getElementProperty(atypes, myATID, 'c_ident', me)
170 >    
171 >    myATID = getFirstMatchingElement(atypes, 'c_ident', c_ident)
172 > !    write(*,*) 'myATID= ', myATID, ' c_ident = ', c_ident
173  
174 <    ShapeMap%atidToShape(me)                         = current
175 <    ShapeMap%Shapes(current)%atid                    = me
174 >    ShapeMap%atidToShape(myATID)                     = current
175 >    ShapeMap%Shapes(current)%atid                    = myATID
176      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
177      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
178      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 336 | Line 338 | contains  
338      integer :: status
339      integer :: i, j, l, m, lm, function_type
340      real(kind=dp) :: thisDP, sigma
341 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
341 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
342      logical :: thisProperty
343  
344      status = 0
# Line 351 | Line 353 | contains  
353      if (nAtypes == 0) then
354         status = -1
355         return
356 <    end if
356 >    end if      
357  
358      ! atypes comes from c side
359 <    do i = 0, nAtypes
360 <
361 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
362 <
359 >    do i = 1, nAtypes
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 >      
365         if (thisProperty) then
366 <
366 >                
367            ShapeMap%currentShape = ShapeMap%currentShape + 1
368            current = ShapeMap%currentShape
369  
370 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
371 <          ShapeMap%atidToShape(thisIP) = current
368 <          ShapeMap%Shapes(current)%atid = thisIP
370 >          ShapeMap%atidToShape(myATID) = current
371 >          ShapeMap%Shapes(current)%atid = myATID
372  
373            ShapeMap%Shapes(current)%isLJ = .true.
374  
375 <          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
376 <          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
375 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
376 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
377  
378         endif
379  
# Line 378 | Line 381 | contains  
381  
382      haveShapeMap = .true.
383  
384 + !    do i = 1, ShapeMap%n_shapes
385 + !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
386 + !    end do
387 +
388    end subroutine complete_Shape_FF
389  
390    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
# Line 492 | Line 499 | contains  
499  
500      !! We assume that the rotation matrices have already been calculated
501      !! and placed in the A array.
495 –
502      r3 = r2*rij
503      r5 = r3*r2
504  
# Line 516 | Line 522 | contains  
522      ! use the atid to find the shape type (st) for each atom:
523      st1 = ShapeMap%atidToShape(atid1)
524      st2 = ShapeMap%atidToShape(atid2)
525 +    
526 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
527  
528      if (ShapeMap%Shapes(st1)%isLJ) then
529  
# Line 666 | Line 674 | contains  
674            endif
675  
676            sigma_i = sigma_i + plm_i(m,l)*Phunc
677 <
677 >          write(*,*) 'dsigmaidux = ', dsigmaidux
678 >          write(*,*) 'Phunc = ', Phunc
679            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
680                 Phunc * dlm_i(m,l) * dctidx
681            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
682                 Phunc * dlm_i(m,l) * dctidy
683            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
684                 Phunc * dlm_i(m,l) * dctidz
676 –
685            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
686                 Phunc * dlm_i(m,l) * dctidux
687            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
688                 Phunc * dlm_i(m,l) * dctiduy
689            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
690                 Phunc * dlm_i(m,l) * dctiduz
691 <
691 >          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
692 >                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
693 >                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
694         end do
695  
696         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 864 | Line 874 | contains  
874         dspjduz = 0.0d0
875  
876  
877 <       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
878 <       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
877 > !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
878 > !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
879         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
880              ShapeMap%Shapes(st2)%bigL, LMAX, &
881              plm_j, dlm_j)
# Line 1005 | Line 1015 | contains  
1015               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1016            endif
1017  
1018 <          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1018 > !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1019  
1020            eps_j = eps_j + plm_j(m,l)*Phunc
1021  
# Line 1030 | Line 1040 | contains  
1040      ! phew, now let's assemble the potential energy:
1041  
1042      sigma = 0.5*(sigma_i + sigma_j)
1043 <
1043 > !    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1044      dsigmadxi = 0.5*dsigmaidx
1045      dsigmadyi = 0.5*dsigmaidy
1046      dsigmadzi = 0.5*dsigmaidz
# Line 1062 | Line 1072 | contains  
1072      dsduzj = 0.5*dsjduz
1073  
1074      eps = sqrt(eps_i * eps_j)
1075 <
1075 >    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1076 >    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1077 > !    write(*,*) sigma_i, ' is sigma i; ', s_i, ' is s i; ', eps_i, ' is eps i'
1078      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1079      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1080      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1077 | Line 1089 | contains  
1089      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1090      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1091  
1092 < !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 < !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1092 > !    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 > !    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1094 >
1095 > !    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1096 > !    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1097   !!$
1083 – !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1084 – !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1085 – !!$
1098   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1099  
1100      rtdenom = rij-sigma+s
# Line 1092 | Line 1104 | contains  
1104      drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1105      drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1106      drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1107 +    write(*,*) dsduxi, ' is dsduxi; ', drduxi, ' is drduxi; ', dsigmaduxi, &
1108 +               ' is dsigmaduxi; ', dsduxi, ' is dsduxi'
1109      drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1110      drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1111      drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
# Line 1101 | Line 1115 | contains  
1115      drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1116      drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1117  
1118 + !    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1119 + !    write(*,*) 'drtdu_i = ', drtduxi, drtduyi, drtduzi
1120 +
1121      rt2 = rt*rt
1122      rt3 = rt2*rt
1123      rt5 = rt2*rt3
# Line 1136 | Line 1153 | contains  
1153      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1154      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1155      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1156 <
1156 >    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1157      ! do the torques first since they are easy:
1158      ! remember that these are still in the body fixed axes
1159  
# Line 1145 | Line 1162 | contains  
1162   !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1163   !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1164   !!$
1165 <    txi =  (dvduzi - dvduyi) * sw
1166 <    tyi =  (dvduxi - dvduzi) * sw
1167 <    tzi =  (dvduyi - dvduxi) * sw
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
1169 > !    txj = (dvduzj - dvduyj) * sw
1170 > !    tyj = (dvduxj - dvduzj) * sw
1171 > !    tzj = (dvduyj - dvduxj) * sw
1172  
1173 < !!$    txi = -dvduxi * sw
1174 < !!$    tyi = -dvduyi * sw
1175 < !!$    tzi = -dvduzi * sw
1176 < !!$
1177 < !!$    txj = dvduxj * sw
1178 < !!$    tyj = dvduyj * sw
1179 < !!$    tzj = dvduzj * sw
1173 >    txi = dvduxi * sw
1174 >    tyi = dvduyi * sw
1175 >    tzi = dvduzi * sw
1176 >
1177 >    txj = dvduxj * sw
1178 >    tyj = dvduyj * sw
1179 >    tzj = dvduzj * sw
1180  
1181      write(*,*) 't1 = ', txi, tyi, tzi
1182      write(*,*) 't2 = ', txj, tyj, tzj
# Line 1201 | Line 1218 | contains  
1218      fyj = dvdyj * sw
1219      fzj = dvdzj * sw
1220  
1221 +
1222   #ifdef IS_MPI
1223      fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1224      fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
# Line 1230 | Line 1248 | contains  
1248      fxradial = 0.5_dp * (fxii + fxji)
1249      fyradial = 0.5_dp * (fyii + fyji)
1250      fzradial = 0.5_dp * (fzii + fzji)
1251 <
1251 >    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1252   #ifdef IS_MPI
1253      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1254      f_Row(2,atom1) = f_Row(2,atom1) + fyradial

Diff Legend

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