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 2323 by chuckv, Fri Sep 23 20:31:02 2005 UTC vs.
Revision 2355 by chuckv, Wed Oct 12 18:59:16 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.49 2005-09-23 20:31:02 chuckv Exp $, $Date: 2005-09-23 20:31:02 $, $Name: not supported by cvs2svn $, $Revision: 1.49 $
48 > !! @version $Id: doForces.F90,v 1.54 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
49  
50  
51   module doForces
# 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 180 | 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 260 | 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 276 | 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 332 | 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 365 | 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 424 | 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 630 | Line 749 | contains
749  
750      !! Stress Tensor
751      real( kind = dp), dimension(9) :: tau  
752 <    real ( kind = dp ) :: pot
752 >    real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
753      logical ( kind = 2) :: do_pot_c, do_stress_c
754      logical :: do_pot
755      logical :: do_stress
756      logical :: in_switching_region
757   #ifdef IS_MPI
758 <    real( kind = DP ) :: pot_local
758 >    real( kind = DP ), dimension(POT_ARRAY_SIZE) :: pot_local
759      integer :: nAtomsInRow
760      integer :: nAtomsInCol
761      integer :: nprocs
# Line 788 | Line 907 | contains
907                    q_group(:,j), d_grp, rgrpsq)
908   #endif      
909  
910 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 909 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 1019 | Line 1139 | contains
1139                  !! potential and torques:
1140                  call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1141   #ifdef IS_MPI
1142 <                pot_local = pot_local + rfpot
1142 >                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1143   #else
1144 <                pot = pot + rfpot
1144 >                pot(RF_POT) = pot(RF_POT) + rfpot
1145  
1146   #endif
1147               endif
# Line 1033 | Line 1153 | contains
1153   #ifdef IS_MPI
1154  
1155      if (do_pot) then
1156 <       pot = pot + pot_local
1156 >       pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1157 >            + pot_local(1:SIZE_POT_ARRAY)
1158         !! we assume the c code will do the allreduce to get the total potential
1159         !! we could do it right here if we needed to...
1160      endif
# Line 1059 | Line 1180 | contains
1180    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1181         eFrame, A, f, t, pot, vpair, fpair)
1182  
1183 <    real( kind = dp ) :: pot, vpair, sw
1183 >    real( kind = dp ) :: vpair, sw
1184 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1185      real( kind = dp ), dimension(3) :: fpair
1186      real( kind = dp ), dimension(nLocal)   :: mfact
1187      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1091 | Line 1213 | contains
1213      iHash = InteractionHash(me_i, me_j)
1214  
1215      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1216 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1216 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1217      endif
1218  
1219      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1220         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1221 <            pot, eFrame, f, t, do_pot)
1221 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1222  
1223         if (electrostaticSummationMethod == REACTION_FIELD) then
1224  
# Line 1109 | Line 1231 | contains
1231  
1232      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1233         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1234 <            pot, A, f, t, do_pot)
1234 >            pot(STICKY_POT), A, f, t, do_pot)
1235      endif
1236  
1237      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1238         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 <            pot, A, f, t, do_pot)
1239 >            pot(STICKYPOWER_POT), A, f, t, do_pot)
1240      endif
1241  
1242      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1243         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1244 <            pot, A, f, t, do_pot)
1244 >            pot(GAYBERNE_POT), A, f, t, do_pot)
1245      endif
1246      
1247      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1248   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1249 < !           pot, A, f, t, do_pot)
1249 > !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1250      endif
1251  
1252      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1253 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1253 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1254              do_pot)
1255      endif
1256  
1257      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1258         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1259 <            pot, A, f, t, do_pot)
1259 >            pot(SHAPE_POT), A, f, t, do_pot)
1260      endif
1261  
1262      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1263         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1264 <            pot, A, f, t, do_pot)
1264 >            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1265      endif
1266      
1267    end subroutine do_pair
# Line 1147 | Line 1269 | contains
1269    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1270         do_pot, do_stress, eFrame, A, f, t, pot)
1271  
1272 <    real( kind = dp ) :: pot, sw
1272 >    real( kind = dp ) :: sw
1273 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1274      real( kind = dp ), dimension(9,nLocal) :: eFrame
1275      real (kind=dp), dimension(9,nLocal) :: A
1276      real (kind=dp), dimension(3,nLocal) :: f
# Line 1182 | Line 1305 | contains
1305  
1306    subroutine do_preforce(nlocal,pot)
1307      integer :: nlocal
1308 <    real( kind = dp ) :: pot
1308 >    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1309  
1310      if (FF_uses_EAM .and. SIM_uses_EAM) then
1311 <       call calc_EAM_preforce_Frho(nlocal,pot)
1311 >       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1312      endif
1313  
1314  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines