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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2317 by chrisfen, Wed Sep 21 17:20:14 2005 UTC vs.
Revision 2375 by gezelter, Mon Oct 17 19:12:45 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.48 2005-09-21 17:20:10 chrisfen Exp $, $Date: 2005-09-21 17:20:10 $, $Name: not supported by cvs2svn $, $Revision: 1.48 $
48 > !! @version $Id: doForces.F90,v 1.59 2005-10-17 19:12:34 gezelter Exp $, $Date: 2005-10-17 19:12:34 $, $Name: not supported by cvs2svn $, $Revision: 1.59 $
49  
50  
51   module doForces
# Line 59 | Line 59 | module doForces
59    use sticky
60    use electrostatic_module
61    use reaction_field_module
62 <  use gb_pair
62 >  use gayberne
63    use shapes
64    use vector_class
65    use eam
# Line 124 | Line 124 | module doForces
124    ! Bit hash to determine pair-pair interactions.
125    integer, dimension(:,:), allocatable :: InteractionHash
126    real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 <  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
128 <  integer, dimension(:), allocatable :: groupToGtype
129 <  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
127 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129 >
130 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
132 >
133 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135    type ::gtypeCutoffs
136       real(kind=dp) :: rcut
137       real(kind=dp) :: rcutsq
# Line 136 | Line 141 | module doForces
141  
142    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143    real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 +  real(kind=dp),save :: listSkin
145    
146   contains
147  
# Line 179 | Line 185 | contains
185  
186      if (.not. allocated(InteractionHash)) then
187         allocate(InteractionHash(nAtypes,nAtypes))
188 +    else
189 +       deallocate(InteractionHash)
190 +       allocate(InteractionHash(nAtypes,nAtypes))
191      endif
192  
193      if (.not. allocated(atypeMaxCutoff)) then
194         allocate(atypeMaxCutoff(nAtypes))
195 +    else
196 +       deallocate(atypeMaxCutoff)
197 +       allocate(atypeMaxCutoff(nAtypes))
198      endif
199          
200      do i = 1, nAtypes
# Line 259 | Line 271 | contains
271      logical :: GtypeFound
272  
273      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
274 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
274 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275      integer :: nGroupsInRow
276 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
276 >    integer :: nGroupsInCol
277 >    integer :: nGroupTypesRow,nGroupTypesCol
278 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279      real(kind=dp) :: biggestAtypeCutoff
280  
281      stat = 0
# Line 275 | Line 289 | contains
289      endif
290   #ifdef IS_MPI
291      nGroupsInRow = getNgroupsInRow(plan_group_row)
292 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
293   #endif
294      nAtypes = getSize(atypes)
295   ! Set all of the initial cutoffs to zero.
# Line 331 | Line 346 | contains
346         endif
347      enddo
348    
349 <    nGroupTypes = 0
349 >
350      
351      istart = 1
352 +    jstart = 1
353   #ifdef IS_MPI
354      iend = nGroupsInRow
355 +    jend = nGroupsInCol
356   #else
357      iend = nGroups
358 +    jend = nGroups
359   #endif
360      
361      !! allocate the groupToGtype and gtypeMaxCutoff here.
362 <    if(.not.allocated(groupToGtype)) then
363 <       allocate(groupToGtype(iend))
364 <       allocate(groupMaxCutoff(iend))
365 <       allocate(gtypeMaxCutoff(iend))
366 <       groupMaxCutoff = 0.0_dp
367 <       gtypeMaxCutoff = 0.0_dp
362 >    if(.not.allocated(groupToGtypeRow)) then
363 >     !  allocate(groupToGtype(iend))
364 >       allocate(groupToGtypeRow(iend))
365 >    else
366 >       deallocate(groupToGtypeRow)
367 >       allocate(groupToGtypeRow(iend))
368      endif
369 +    if(.not.allocated(groupMaxCutoffRow)) then
370 +       allocate(groupMaxCutoffRow(iend))
371 +    else
372 +       deallocate(groupMaxCutoffRow)
373 +       allocate(groupMaxCutoffRow(iend))
374 +    end if
375 +
376 +    if(.not.allocated(gtypeMaxCutoffRow)) then
377 +       allocate(gtypeMaxCutoffRow(iend))
378 +    else
379 +       deallocate(gtypeMaxCutoffRow)
380 +       allocate(gtypeMaxCutoffRow(iend))
381 +    endif
382 +
383 +
384 + #ifdef IS_MPI
385 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
386 +    if(.not.associated(groupToGtypeCol)) then
387 +       allocate(groupToGtypeCol(jend))
388 +    else
389 +       deallocate(groupToGtypeCol)
390 +       allocate(groupToGtypeCol(jend))
391 +    end if
392 +
393 +    if(.not.associated(groupToGtypeCol)) then
394 +       allocate(groupToGtypeCol(jend))
395 +    else
396 +       deallocate(groupToGtypeCol)
397 +       allocate(groupToGtypeCol(jend))
398 +    end if
399 +    if(.not.associated(gtypeMaxCutoffCol)) then
400 +       allocate(gtypeMaxCutoffCol(jend))
401 +    else
402 +       deallocate(gtypeMaxCutoffCol)      
403 +       allocate(gtypeMaxCutoffCol(jend))
404 +    end if
405 +
406 +       groupMaxCutoffCol = 0.0_dp
407 +       gtypeMaxCutoffCol = 0.0_dp
408 +
409 + #endif
410 +       groupMaxCutoffRow = 0.0_dp
411 +       gtypeMaxCutoffRow = 0.0_dp
412 +
413 +
414      !! first we do a single loop over the cutoff groups to find the
415      !! largest cutoff for any atypes present in this group.  We also
416      !! create gtypes at this point.
417      
418      tol = 1.0d-6
419 <    
419 >    nGroupTypesRow = 0
420 >
421      do i = istart, iend      
422         n_in_i = groupStartRow(i+1) - groupStartRow(i)
423 <       groupMaxCutoff(i) = 0.0_dp
423 >       groupMaxCutoffRow(i) = 0.0_dp
424         do ia = groupStartRow(i), groupStartRow(i+1)-1
425            atom1 = groupListRow(ia)
426   #ifdef IS_MPI
# Line 364 | Line 428 | contains
428   #else
429            me_i = atid(atom1)
430   #endif          
431 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
432 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
431 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433            endif          
434         enddo
435  
436 <       if (nGroupTypes.eq.0) then
437 <          nGroupTypes = nGroupTypes + 1
438 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
439 <          groupToGtype(i) = nGroupTypes
436 >       if (nGroupTypesRow.eq.0) then
437 >          nGroupTypesRow = nGroupTypesRow + 1
438 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439 >          groupToGtypeRow(i) = nGroupTypesRow
440         else
441            GtypeFound = .false.
442 <          do g = 1, nGroupTypes
443 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
444 <                groupToGtype(i) = g
442 >          do g = 1, nGroupTypesRow
443 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444 >                groupToGtypeRow(i) = g
445                  GtypeFound = .true.
446               endif
447            enddo
448            if (.not.GtypeFound) then            
449 <             nGroupTypes = nGroupTypes + 1
450 <             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
451 <             groupToGtype(i) = nGroupTypes
449 >             nGroupTypesRow = nGroupTypesRow + 1
450 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451 >             groupToGtypeRow(i) = nGroupTypesRow
452            endif
453         endif
454      enddo    
455  
456 + #ifdef IS_MPI
457 +    do j = jstart, jend      
458 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
459 +       groupMaxCutoffCol(j) = 0.0_dp
460 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
461 +          atom1 = groupListCol(ja)
462 +
463 +          me_j = atid_col(atom1)
464 +
465 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467 +          endif          
468 +       enddo
469 +
470 +       if (nGroupTypesCol.eq.0) then
471 +          nGroupTypesCol = nGroupTypesCol + 1
472 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473 +          groupToGtypeCol(j) = nGroupTypesCol
474 +       else
475 +          GtypeFound = .false.
476 +          do g = 1, nGroupTypesCol
477 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478 +                groupToGtypeCol(j) = g
479 +                GtypeFound = .true.
480 +             endif
481 +          enddo
482 +          if (.not.GtypeFound) then            
483 +             nGroupTypesCol = nGroupTypesCol + 1
484 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485 +             groupToGtypeCol(j) = nGroupTypesCol
486 +          endif
487 +       endif
488 +    enddo    
489 +
490 + #else
491 + ! Set pointers to information we just found
492 +    nGroupTypesCol = nGroupTypesRow
493 +    groupToGtypeCol => groupToGtypeRow
494 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
495 +    groupMaxCutoffCol => groupMaxCutoffRow
496 + #endif
497 +
498 +
499 +
500 +
501 +
502      !! allocate the gtypeCutoffMap here.
503 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
503 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504      !! then we do a double loop over all the group TYPES to find the cutoff
505      !! map between groups of two types
506 <    
507 <    do i = 1, nGroupTypes
508 <       do j = 1, nGroupTypes
506 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 >
508 >    do i = 1, nGroupTypesRow
509 >       do j = 1, nGroupTypesCol
510        
511            select case(cutoffPolicy)
512            case(TRADITIONAL_CUTOFF_POLICY)
513 <             thisRcut = maxval(gtypeMaxCutoff)
513 >             thisRcut = tradRcut
514            case(MIX_CUTOFF_POLICY)
515 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
515 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516            case(MAX_CUTOFF_POLICY)
517 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
517 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518            case default
519               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520               return
# Line 411 | Line 522 | contains
522            gtypeCutoffMap(i,j)%rcut = thisRcut
523            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524            skin = defaultRlist - defaultRcut
525 +          listSkin = skin ! set neighbor list skin thickness
526            gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527  
528            ! sanity check
# Line 422 | Line 534 | contains
534            endif
535         enddo
536      enddo
537 <
537 >    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538 >    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539 >    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540 > #ifdef IS_MPI
541 >    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542 >    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543 > #endif
544 >    groupMaxCutoffCol => null()
545 >    gtypeMaxCutoffCol => null()
546 >    
547      haveGtypeCutoffMap = .true.
548     end subroutine createGtypeCutoffMap
549  
# Line 586 | Line 707 | contains
707         end if
708      endif
709  
589    if (FF_uses_GayBerne) then
590       call check_gb_pair_FF(my_status)
591       if (my_status .ne. 0) then
592          thisStat = -1
593          haveSaneForceField = .false.
594          return
595       endif
596    endif
597
710      if (.not. haveNeighborList) then
711         !! Create neighbor lists
712         call expandNeighborList(nLocal, my_status)
# Line 628 | Line 740 | contains
740  
741      !! Stress Tensor
742      real( kind = dp), dimension(9) :: tau  
743 <    real ( kind = dp ) :: pot
743 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
744      logical ( kind = 2) :: do_pot_c, do_stress_c
745      logical :: do_pot
746      logical :: do_stress
747      logical :: in_switching_region
748   #ifdef IS_MPI
749 <    real( kind = DP ) :: pot_local
749 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
750      integer :: nAtomsInRow
751      integer :: nAtomsInCol
752      integer :: nprocs
# Line 659 | Line 771 | contains
771      integer :: propPack_i, propPack_j
772      integer :: loopStart, loopEnd, loop
773      integer :: iHash
774 <    real(kind=dp) :: listSkin = 1.0  
774 >  
775  
776      !! initialize local variables  
777  
# Line 786 | Line 898 | contains
898                    q_group(:,j), d_grp, rgrpsq)
899   #endif      
900  
901 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
901 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
902                  if (update_nlist) then
903                     nlist = nlist + 1
904  
# Line 907 | Line 1019 | contains
1019                  endif
1020               end if
1021            enddo
1022 +
1023         enddo outer
1024  
1025         if (update_nlist) then
# Line 966 | Line 1079 | contains
1079  
1080      if (do_pot) then
1081         ! scatter/gather pot_row into the members of my column
1082 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1083 <
1082 >       do i = 1,LR_POT_TYPES
1083 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1084 >       end do
1085         ! scatter/gather pot_local into all other procs
1086         ! add resultant to get total pot
1087         do i = 1, nlocal
1088 <          pot_local = pot_local + pot_Temp(i)
1088 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1089 >               + pot_Temp(1:LR_POT_TYPES,i)
1090         enddo
1091  
1092         pot_Temp = 0.0_DP
1093 <
1094 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1093 >       do i = 1,LR_POT_TYPES
1094 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1095 >       end do
1096         do i = 1, nlocal
1097 <          pot_local = pot_local + pot_Temp(i)
1097 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1098 >               + pot_Temp(1:LR_POT_TYPES,i)
1099         enddo
1100  
1101      endif
# Line 1017 | Line 1134 | contains
1134                  !! potential and torques:
1135                  call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1136   #ifdef IS_MPI
1137 <                pot_local = pot_local + rfpot
1137 >                pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1138   #else
1139 <                pot = pot + rfpot
1139 >                pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1140  
1141   #endif
1142               endif
# Line 1031 | Line 1148 | contains
1148   #ifdef IS_MPI
1149  
1150      if (do_pot) then
1151 <       pot = pot + pot_local
1151 >       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1152 >            + pot_local(1:LR_POT_TYPES)
1153         !! we assume the c code will do the allreduce to get the total potential
1154         !! we could do it right here if we needed to...
1155      endif
# Line 1057 | Line 1175 | contains
1175    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1176         eFrame, A, f, t, pot, vpair, fpair)
1177  
1178 <    real( kind = dp ) :: pot, vpair, sw
1178 >    real( kind = dp ) :: vpair, sw
1179 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1180      real( kind = dp ), dimension(3) :: fpair
1181      real( kind = dp ), dimension(nLocal)   :: mfact
1182      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1089 | Line 1208 | contains
1208      iHash = InteractionHash(me_i, me_j)
1209  
1210      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1211 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1211 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1212 >            pot(VDW_POT), f, do_pot)
1213      endif
1214  
1215      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1216         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1217 <            pot, eFrame, f, t, do_pot)
1217 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1218  
1219         if (electrostaticSummationMethod == REACTION_FIELD) then
1220  
# Line 1107 | Line 1227 | contains
1227  
1228      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1229         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230 <            pot, A, f, t, do_pot)
1230 >            pot(HB_POT), A, f, t, do_pot)
1231      endif
1232  
1233      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1234         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235 <            pot, A, f, t, do_pot)
1235 >            pot(HB_POT), A, f, t, do_pot)
1236      endif
1237  
1238      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1239         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1240 <            pot, A, f, t, do_pot)
1240 >            pot(VDW_POT), A, f, t, do_pot)
1241      endif
1242      
1243      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1244 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245 < !           pot, A, f, t, do_pot)
1244 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245 >            pot(VDW_POT), A, f, t, do_pot)
1246      endif
1247  
1248      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1249 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1250 <            do_pot)
1249 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1250 >            pot(METALLIC_POT), f, do_pot)
1251      endif
1252  
1253      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1254         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1255 <            pot, A, f, t, do_pot)
1255 >            pot(VDW_POT), A, f, t, do_pot)
1256      endif
1257  
1258      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1259         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1260 <            pot, A, f, t, do_pot)
1260 >            pot(VDW_POT), A, f, t, do_pot)
1261      endif
1262      
1263    end subroutine do_pair
# Line 1145 | Line 1265 | contains
1265    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1266         do_pot, do_stress, eFrame, A, f, t, pot)
1267  
1268 <    real( kind = dp ) :: pot, sw
1268 >    real( kind = dp ) :: sw
1269 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1270      real( kind = dp ), dimension(9,nLocal) :: eFrame
1271      real (kind=dp), dimension(9,nLocal) :: A
1272      real (kind=dp), dimension(3,nLocal) :: f
# Line 1180 | Line 1301 | contains
1301  
1302    subroutine do_preforce(nlocal,pot)
1303      integer :: nlocal
1304 <    real( kind = dp ) :: pot
1304 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1305  
1306      if (FF_uses_EAM .and. SIM_uses_EAM) then
1307 <       call calc_EAM_preforce_Frho(nlocal,pot)
1307 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1308      endif
1309  
1310  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines