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 2432 by chuckv, Tue Nov 15 16:01:06 2005 UTC vs.
Revision 2461 by gezelter, Mon Nov 21 22:59:02 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.68 2005-11-15 16:01:06 chuckv Exp $, $Date: 2005-11-15 16:01:06 $, $Name: not supported by cvs2svn $, $Revision: 1.68 $
48 > !! @version $Id: doForces.F90,v 1.69 2005-11-21 22:58:35 gezelter Exp $, $Date: 2005-11-21 22:58:35 $, $Name: not supported by cvs2svn $, $Revision: 1.69 $
49  
50  
51   module doForces
# Line 87 | Line 87 | module doForces
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89    logical, save :: haveDefaultCutoffs = .false.
90 <  logical, save :: haveRlist = .false.
90 >  logical, save :: haveSkinThickness = .false.
91 >  logical, save :: haveElectrostaticSummationMethod = .false.
92 >  logical, save :: haveCutoffPolicy = .false.
93 >  logical, save :: VisitCutoffsAfterComputing = .false.
94  
95    logical, save :: FF_uses_DirectionalAtoms
96    logical, save :: FF_uses_Dipoles
# Line 106 | Line 109 | module doForces
109    logical, save :: SIM_uses_PBC
110  
111    integer, save :: electrostaticSummationMethod
112 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
113  
114 +  real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
115 +  real(kind=dp), save :: skinThickness
116 +  logical, save :: defaultDoShift
117 +
118    public :: init_FF
119 <  public :: setDefaultCutoffs
119 >  public :: setCutoffs
120 >  public :: cWasLame
121 >  public :: setElectrostaticMethod
122 >  public :: setCutoffPolicy
123 >  public :: setSkinThickness
124    public :: do_force_loop
113  public :: createInteractionHash
114  public :: createGtypeCutoffMap
115  public :: getStickyCut
116  public :: getStickyPowerCut
117  public :: getGayBerneCut
118  public :: getEAMCut
119  public :: getShapeCut
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 144 | Line 149 | module doForces
149    end type gtypeCutoffs
150    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151  
147  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
148  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
149  real(kind=dp),save :: listSkin
150  
152   contains
153  
154 <  subroutine createInteractionHash(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
155    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158      integer :: iHash
# Line 177 | Line 177 | contains
177      logical :: j_is_MEAM
178      real(kind=dp) :: myRcut
179  
180
181    status = 0  
182
180      if (.not. associated(atypes)) then
181 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
185 <       status = -1
181 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182         return
183      endif
184      
185      nAtypes = getSize(atypes)
186      
187      if (nAtypes == 0) then
188 <       status = -1
188 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
189         return
190      end if
191  
# Line 276 | Line 272 | contains
272      haveInteractionHash = .true.
273    end subroutine createInteractionHash
274  
275 <  subroutine createGtypeCutoffMap(stat)
275 >  subroutine createGtypeCutoffMap()
276  
281    integer, intent(out), optional :: stat
277      logical :: i_is_LJ
278      logical :: i_is_Elect
279      logical :: i_is_Sticky
# Line 293 | Line 288 | contains
288      integer :: nGroupsInRow
289      integer :: nGroupsInCol
290      integer :: nGroupTypesRow,nGroupTypesCol
291 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
291 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
292      real(kind=dp) :: biggestAtypeCutoff
293  
299    stat = 0
294      if (.not. haveInteractionHash) then
295 <       call createInteractionHash(myStatus)      
302 <       if (myStatus .ne. 0) then
303 <          write(default_error, *) 'createInteractionHash failed in doForces!'
304 <          stat = -1
305 <          return
306 <       endif
295 >       call createInteractionHash()      
296      endif
297   #ifdef IS_MPI
298      nGroupsInRow = getNgroupsInRow(plan_group_row)
# Line 355 | Line 344 | contains
344                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345               endif
346            endif
347 <          
359 <          
347 >                    
348            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349               biggestAtypeCutoff = atypeMaxCutoff(i)
350            endif
351  
352         endif
353      enddo
366  
367
354      
355      istart = 1
356      jstart = 1
# Line 513 | Line 499 | contains
499      groupMaxCutoffCol => groupMaxCutoffRow
500   #endif
501  
516
517
518
519
502      !! allocate the gtypeCutoffMap here.
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      tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507  
508 <    do i = 1, nGroupTypesRow
508 >    do i = 1, nGroupTypesRow      
509         do j = 1, nGroupTypesCol
510        
511            select case(cutoffPolicy)
# Line 538 | Line 520 | contains
520               return
521            end select
522            gtypeCutoffMap(i,j)%rcut = thisRcut
523 +          
524 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
525 +
526            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
542          skin = defaultRlist - defaultRcut
543          listSkin = skin ! set neighbor list skin thickness
544          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527  
528 +          if (.not.haveSkinThickness) then
529 +             skinThickness = 1.0_dp
530 +          endif
531 +
532 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
533 +
534            ! sanity check
535  
536            if (haveDefaultCutoffs) then
# Line 552 | Line 540 | contains
540            endif
541         enddo
542      enddo
543 +
544      if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
545      if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
546      if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
# Line 565 | Line 554 | contains
554      haveGtypeCutoffMap = .true.
555     end subroutine createGtypeCutoffMap
556  
557 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
569 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
570 <     integer, intent(in) :: cutPolicy
557 >   subroutine setCutoffs(defRcut, defRsw)
558  
559 +     real(kind=dp),intent(in) :: defRcut, defRsw
560 +     character(len = statusMsgSize) :: errMsg
561 +     integer :: localError
562 +
563       defaultRcut = defRcut
564       defaultRsw = defRsw
565 <     defaultRlist = defRlist
566 <     cutoffPolicy = cutPolicy
565 >    
566 >     defaultDoShift = .false.
567 >     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
568 >        
569 >        write(errMsg, *) &
570 >             'cutoffRadius and switchingRadius are set to the same', newline &
571 >             // tab, 'value.  OOPSE will use shifted ', newline &
572 >             // tab, 'potentials instead of switching functions.'
573 >        
574 >        call handleInfo("setCutoffs", errMsg)
575 >        
576 >        defaultDoShift = .true.
577 >        
578 >     endif
579  
580 +     localError = 0
581 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
582 +     call setCutoffEAM( defaultRcut, localError)
583 +     if (localError /= 0) then
584 +       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
585 +       call handleError("setCutoffs", errMsg)
586 +     end if
587 +     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
588 +    
589       haveDefaultCutoffs = .true.
590 <   end subroutine setDefaultCutoffs
590 >   end subroutine setCutoffs
591  
592 +   subroutine cWasLame()
593 +    
594 +     VisitCutoffsAfterComputing = .true.
595 +     return
596 +    
597 +   end subroutine cWasLame
598 +  
599     subroutine setCutoffPolicy(cutPolicy)
600 <
600 >    
601       integer, intent(in) :: cutPolicy
602 +    
603       cutoffPolicy = cutPolicy
604 +     haveCutoffPolicy = .true.
605 +
606       call createGtypeCutoffMap()
585   end subroutine setCutoffPolicy
586    
607      
608 <  subroutine setSimVariables()
609 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
610 <    SIM_uses_EAM = SimUsesEAM()
591 <    SIM_uses_SC  = SimUsesSC()
592 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
593 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
594 <    SIM_uses_PBC = SimUsesPBC()
608 >   end subroutine setCutoffPolicy
609 >  
610 >   subroutine setElectrostaticMethod( thisESM )
611  
612 <    haveSIMvariables = .true.
612 >     integer, intent(in) :: thisESM
613  
614 <    return
615 <  end subroutine setSimVariables
614 >     electrostaticSummationMethod = thisESM
615 >     haveElectrostaticSummationMethod = .true.
616 >    
617 >   end subroutine setElectrostaticMethod
618  
619 +   subroutine setSkinThickness( thisSkin )
620 +    
621 +     real(kind=dp), intent(in) :: thisSkin
622 +    
623 +     skinThickness = thisSkin
624 +     haveSkinThickness = .true.
625 +    
626 +     call createGtypeCutoffMap()
627 +    
628 +   end subroutine setSkinThickness
629 +      
630 +   subroutine setSimVariables()
631 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
632 +     SIM_uses_EAM = SimUsesEAM()
633 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
634 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
635 +     SIM_uses_PBC = SimUsesPBC()
636 +    
637 +     haveSIMvariables = .true.
638 +    
639 +     return
640 +   end subroutine setSimVariables
641 +
642    subroutine doReadyCheck(error)
643      integer, intent(out) :: error
644  
# Line 606 | Line 647 | contains
647      error = 0
648  
649      if (.not. haveInteractionHash) then      
650 <       myStatus = 0      
610 <       call createInteractionHash(myStatus)      
611 <       if (myStatus .ne. 0) then
612 <          write(default_error, *) 'createInteractionHash failed in doForces!'
613 <          error = -1
614 <          return
615 <       endif
650 >       call createInteractionHash()      
651      endif
652  
653      if (.not. haveGtypeCutoffMap) then        
654 <       myStatus = 0      
620 <       call createGtypeCutoffMap(myStatus)      
621 <       if (myStatus .ne. 0) then
622 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
623 <          error = -1
624 <          return
625 <       endif
654 >       call createGtypeCutoffMap()      
655      endif
656  
657 +
658 +    if (VisitCutoffsAfterComputing) then
659 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
660 +    endif
661 +
662 +
663      if (.not. haveSIMvariables) then
664         call setSimVariables()
665      endif
# Line 658 | Line 693 | contains
693    end subroutine doReadyCheck
694  
695  
696 <  subroutine init_FF(thisESM, thisStat)
696 >  subroutine init_FF(thisStat)
697  
663    integer, intent(in) :: thisESM
698      integer, intent(out) :: thisStat  
699      integer :: my_status, nMatches
700      integer, pointer :: MatchList(:) => null()
# Line 668 | Line 702 | contains
702      !! assume things are copacetic, unless they aren't
703      thisStat = 0
704  
671    electrostaticSummationMethod = thisESM
672
705      !! init_FF is called *after* all of the atom types have been
706      !! defined in atype_module using the new_atype subroutine.
707      !!
# Line 766 | Line 798 | contains
798      real( kind = DP ) :: rVal
799      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
800      real(kind=dp) :: rfpot, mu_i, virial
801 +    real(kind=dp):: rCut
802      integer :: me_i, me_j, n_in_i, n_in_j
803      logical :: is_dp_i
804      integer :: neighborListSize
# Line 839 | Line 872 | contains
872         ! (but only on the first time through):
873         if (loop .eq. loopStart) then
874   #ifdef IS_MPI
875 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
875 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
876                 update_nlist)
877   #else
878 <          call checkNeighborList(nGroups, q_group, listSkin, &
878 >          call checkNeighborList(nGroups, q_group, skinThickness, &
879                 update_nlist)
880   #endif
881         endif
# Line 922 | Line 955 | contains
955  
956                     list(nlist) = j
957                  endif
958 +
959 +
960                  
961                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
962  
963 +                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
964                     if (loop .eq. PAIR_LOOP) then
965                        vij = 0.0d0
966                        fij(1:3) = 0.0d0
967                     endif
968                    
969 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
970 <                        in_switching_region)
969 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
970 >                        group_switch, in_switching_region)
971                    
972                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
973                    
# Line 961 | Line 997 | contains
997                           if (loop .eq. PREPAIR_LOOP) then
998   #ifdef IS_MPI                      
999                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1000 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1000 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1001                                   eFrame, A, f, t, pot_local)
1002   #else
1003                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1004 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1005                                   eFrame, A, f, t, pot)
1006   #endif                                              
1007                           else
1008   #ifdef IS_MPI                      
1009                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1010                                   do_pot, eFrame, A, f, t, pot_local, vpair, &
1011 <                                 fpair, d_grp, rgrp)
1011 >                                 fpair, d_grp, rgrp, rCut)
1012   #else
1013                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014                                   do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1015 <                                 d_grp, rgrp)
1015 >                                 d_grp, rgrp, rCut)
1016   #endif
1017                              vij = vij + vpair
1018                              fij(1:3) = fij(1:3) + fpair(1:3)
# Line 1184 | Line 1220 | contains
1220    end subroutine do_force_loop
1221  
1222    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1223 <       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1223 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1224  
1225      real( kind = dp ) :: vpair, sw
1226      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
# Line 1201 | Line 1237 | contains
1237      real ( kind = dp ), intent(inout) :: r_grp
1238      real ( kind = dp ), intent(inout) :: d(3)
1239      real ( kind = dp ), intent(inout) :: d_grp(3)
1240 +    real ( kind = dp ), intent(inout) :: rCut
1241      real ( kind = dp ) :: r
1242      integer :: me_i, me_j
1243  
# Line 1221 | Line 1258 | contains
1258      iHash = InteractionHash(me_i, me_j)
1259      
1260      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1261 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1261 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1262              pot(VDW_POT), f, do_pot)
1263      endif
1264      
1265      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1266 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1266 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1267              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1268      endif
1269      
# Line 1246 | Line 1283 | contains
1283      endif
1284      
1285      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1286 <       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1286 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1287              pot(VDW_POT), A, f, t, do_pot)
1288      endif
1289      
# Line 1266 | Line 1303 | contains
1303      endif
1304  
1305      if ( iand(iHash, SC_PAIR).ne.0 ) then      
1306 <       call do_SC_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1306 >       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1307              pot(METALLIC_POT), f, do_pot)
1308      endif
1309  
# Line 1274 | Line 1311 | contains
1311      
1312    end subroutine do_pair
1313  
1314 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1314 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1315         do_pot, do_stress, eFrame, A, f, t, pot)
1316  
1317      real( kind = dp ) :: sw
# Line 1286 | Line 1323 | contains
1323  
1324      logical, intent(inout) :: do_pot, do_stress
1325      integer, intent(in) :: i, j
1326 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1326 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1327      real ( kind = dp )                :: r, rc
1328      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1329  
# Line 1305 | Line 1342 | contains
1342      iHash = InteractionHash(me_i, me_j)
1343  
1344      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1345 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1345 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1346      endif
1347  
1348      if ( iand(iHash, SC_PAIR).ne.0 ) then      
1349 <            call calc_SC_prepair_rho(i, j, d, r, rijsq )
1349 >            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1350      endif
1351      
1352    end subroutine do_prepair

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines