--- trunk/OOPSE/libmdtools/do_Forces.F90 2004/01/05 22:49:14 898 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/01/06 18:54:57 900 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.44 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $ +!! @version $Id: do_Forces.F90,v 1.45 2004-01-06 18:54:57 chuckv Exp $, $Date: 2004-01-06 18:54:57 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $ module do_Forces use force_globals @@ -30,14 +30,28 @@ module do_Forces #define __FORTRAN90 #include "fForceField.h" - logical, save :: do_forces_initialized = .false., haveRlist = .false. + logical, save :: haveRlist = .false. + logical, save :: haveNeighborList = .false. logical, save :: havePolicies = .false. + logical, save :: haveSIMvariables = .false. + logical, save :: havePropertyMap = .false. + logical, save :: haveSaneForceField = .false. logical, save :: FF_uses_LJ logical, save :: FF_uses_sticky logical, save :: FF_uses_dipoles logical, save :: FF_uses_RF logical, save :: FF_uses_GB logical, save :: FF_uses_EAM + logical, save :: SIM_uses_LJ + logical, save :: SIM_uses_sticky + logical, save :: SIM_uses_dipoles + logical, save :: SIM_uses_RF + logical, save :: SIM_uses_GB + logical, save :: SIM_uses_EAM + logical, save :: SIM_requires_postpair_calc + logical, save :: SIM_requires_prepair_calc + logical, save :: SIM_uses_directional_atoms + logical, save :: SIM_uses_PBC real(kind=dp), save :: rlist, rlistsq @@ -45,6 +59,7 @@ module do_Forces public :: do_force_loop public :: setRlistDF + #ifdef PROFILE public :: getforcetime real, save :: forceTime = 0 @@ -52,9 +67,17 @@ module do_Forces integer :: nLoops #endif - logical, allocatable :: propertyMapI(:,:) - logical, allocatable :: propertyMapJ(:,:) + type :: Properties + logical :: is_lj = .false. + logical :: is_sticky = .false. + logical :: is_dp = .false. + logical :: is_gb = .false. + logical :: is_eam = .false. + real(kind=DP) :: dipole_moment = 0.0_DP + end type Properties + type(Properties), dimension(:),allocatable :: PropertyMap + contains subroutine setRlistDF( this_rlist ) @@ -65,10 +88,130 @@ contains rlistsq = rlist * rlist haveRlist = .true. - if( havePolicies ) do_forces_initialized = .true. end subroutine setRlistDF + + subroutine createPropertyMap(status) + integer :: nAtypes + integer :: status + integer :: i + logical :: thisProperty + real (kind=DP) :: thisDPproperty + + status = 0 + + nAtypes = getSize(atypes) + + if (nAtypes == 0) then + status = -1 + return + end if + + if (.not. allocated(PropertyMap)) then + allocate(PropertyMap(nAtypes)) + endif + + do i = 1, nAtypes + call getElementProperty(atypes, i, "is_LJ", thisProperty) + PropertyMap(i)%is_LJ = thisProperty + call getElementProperty(atypes, i, "is_DP", thisProperty) + PropertyMap(i)%is_DP = thisProperty + + if (thisProperty) then + call getElementProperty(atypes, i, "dipole_moment", thisDPproperty) + PropertyMap(i)%dipole_moment = thisDPproperty + endif + + call getElementProperty(atypes, i, "is_Sticky", thisProperty) + PropertyMap(i)%is_Sticky = thisProperty + call getElementProperty(atypes, i, "is_GB", thisProperty) + PropertyMap(i)%is_GB = thisProperty + call getElementProperty(atypes, i, "is_EAM", thisProperty) + PropertyMap(i)%is_EAM = thisProperty + end do + + havePropertyMap = .true. + + end subroutine createPropertyMap + + subroutine setSimVariables() + SIM_uses_LJ = SimUsesLJ() + SIM_uses_sticky = SimUsesSticky() + SIM_uses_dipoles = SimUsesDipoles() + SIM_uses_RF = SimUsesRF() + SIM_uses_GB = SimUsesGB() + SIM_uses_EAM = SimUsesEAM() + SIM_requires_postpair_calc = SimRequiresPostpairCalc() + SIM_requires_prepair_calc = SimRequiresPrepairCalc() + SIM_uses_directional_atoms = SimUsesDirectionalAtoms() + SIM_uses_PBC = SimUsesPBC() + + haveSIMvariables = .true. + return + end subroutine setSimVariables + + subroutine doReadyCheck(error) + integer, intent(out) :: error + + integer :: myStatus + + error = 0 + + if (.not. havePropertyMap) then + + myStatus = 0 + + call createPropertyMap(myStatus) + + if (myStatus .ne. 0) then + write(default_error, *) 'createPropertyMap failed in do_Forces!' + error = -1 + return + endif + endif + + if (.not. haveSIMvariables) then + call setSimVariables() + endif + + if (.not. haveRlist) then + write(default_error, *) 'rList has not been set in do_Forces!' + error = -1 + return + endif + + if (SIM_uses_LJ .and. FF_uses_LJ) then + if (.not. havePolicies) then + write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!' + error = -1 + return + endif + endif + + if (.not. haveNeighborList) then + write(default_error, *) 'neighbor list has not been initialized in do_Forces!' + error = -1 + return + end if + + if (.not. haveSaneForceField) then + write(default_error, *) 'Force Field is not sane in do_Forces!' + error = -1 + return + end if + +#ifdef IS_MPI + if (.not. isMPISimSet()) then + write(default_error,*) "ERROR: mpiSimulation has not been initialized!" + error = -1 + return + endif +#endif + return + end subroutine doReadyCheck + + subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat) integer, intent(in) :: LJMIXPOLICY @@ -113,6 +256,9 @@ contains call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) if (nMatches .gt. 0) FF_uses_EAM = .true. + !! Assume sanity (for the sake of argument) + haveSaneForceField = .true. + !! check to make sure the FF_uses_RF setting makes sense if (FF_uses_dipoles) then @@ -124,6 +270,7 @@ contains if (FF_uses_RF) then write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' thisStat = -1 + haveSaneForceField = .false. return endif endif @@ -138,18 +285,22 @@ contains case default write(default_error,*) 'unknown LJ Mixing Policy!' thisStat = -1 + haveSaneForceField = .false. return end select if (my_status /= 0) then thisStat = -1 + haveSaneForceField = .false. return end if + havePolicies = .true. endif if (FF_uses_sticky) then call check_sticky_FF(my_status) if (my_status /= 0) then thisStat = -1 + haveSaneForceField = .false. return end if endif @@ -158,25 +309,25 @@ contains if (FF_uses_EAM) then call init_EAM_FF(my_status) if (my_status /= 0) then - write(*,*) "init_EAM_FF returned a bad status" + write(default_error, *) "init_EAM_FF returned a bad status" thisStat = -1 + haveSaneForceField = .false. return end if endif - - if (FF_uses_GB) then call check_gb_pair_FF(my_status) if (my_status .ne. 0) then thisStat = -1 + haveSaneForceField = .false. return endif endif if (FF_uses_GB .and. FF_uses_LJ) then endif - if (.not. do_forces_initialized) then + if (.not. haveNeighborList) then !! Create neighbor lists call expandNeighborList(nLocal, my_status) if (my_Status /= 0) then @@ -184,12 +335,9 @@ contains thisStat = -1 return endif + haveNeighborList = .true. endif - - havePolicies = .true. - if( haveRlist ) do_forces_initialized = .true. - end subroutine init_FF @@ -246,9 +394,9 @@ contains natoms = nlocal #endif - call check_initialization(localError) + call doReadyCheck(localError) if ( localError .ne. 0 ) then - call handleError("do_force_loop","Not Initialized") + call handleError("do_force_loop", "Not Initialized") error = -1 return end if @@ -256,78 +404,7 @@ contains do_pot = do_pot_c do_stress = do_stress_c - - -#ifdef IS_MPI - if (.not.allocated(propertyMapI)) then - allocate(propertyMapI(5,nrow)) - endif - - do i = 1, nrow - me_i = atid_row(i) -#else - if (.not.allocated(propertyMapI)) then - allocate(propertyMapI(5,nlocal)) - endif - - do i = 1, natoms - me_i = atid(i) -#endif - - propertyMapI(1:5,i) = .false. - - call getElementProperty(atypes, me_i, "propertyPack", propPack_i) - - ! unpack the properties - - if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) & - propertyMapI(1, i) = .true. - if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) & - propertyMapI(2, i) = .true. - if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) & - propertyMapI(3, i) = .true. - if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) & - propertyMapI(4, i) = .true. - if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) & - propertyMapI(5, i) = .true. - - end do - -#ifdef IS_MPI - if (.not.allocated(propertyMapJ)) then - allocate(propertyMapJ(5,ncol)) - endif - do j = 1, ncol - me_j = atid_col(j) -#else - if (.not.allocated(propertyMapJ)) then - allocate(propertyMapJ(5,nlocal)) - endif - - do j = 1, natoms - me_j = atid(j) -#endif - - propertyMapJ(1:5,j) = .false. - - call getElementProperty(atypes, me_j, "propertyPack", propPack_j) - - ! unpack the properties - - if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) & - propertyMapJ(1, j) = .true. - if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) & - propertyMapJ(2, j) = .true. - if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) & - propertyMapJ(3, j) = .true. - if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) & - propertyMapJ(4, j) = .true. - if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) & - propertyMapJ(5, j) = .true. - - end do - ! Gather all information needed by all force loops: #ifdef IS_MPI @@ -335,7 +412,7 @@ contains call gather(q,q_Row,plan_row3d) call gather(q,q_Col,plan_col3d) - if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then + if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then call gather(u_l,u_l_Row,plan_row3d) call gather(u_l,u_l_Col,plan_col3d) @@ -351,7 +428,7 @@ contains nloops = nloops + 1 #endif - if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then + if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then !! See if we need to update neighbor lists call checkNeighborList(nlocal, q, listSkin, update_nlist) !! if_mpi_gather_stuff_for_prepair @@ -672,7 +749,7 @@ contains f(1:3,i) = f(1:3,i) + f_temp(1:3,i) end do - if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then + if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then t_temp = 0.0_dp call scatter(t_Row,t_temp,plan_row3d) do i = 1,nlocal @@ -706,9 +783,9 @@ contains endif #endif - if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then + if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then - if (FF_uses_RF .and. SimUsesRF()) then + if (FF_uses_RF .and. SIM_uses_RF) then #ifdef IS_MPI call scatter(rf_Row,rf,plan_row3d) @@ -726,12 +803,14 @@ contains #else me_i = atid(i) #endif - call getElementProperty(atypes, me_i, "is_DP", is_DP_i) - if ( is_DP_i ) then - call getElementProperty(atypes, me_i, "dipole_moment", mu_i) + + if (PropertyMap(me_i)%is_DP) then + + mu_i = PropertyMap(me_i)%dipole_moment + !! The reaction field needs to include a self contribution !! to the field: - call accumulate_self_rf(i, mu_i, u_l) + call accumulate_self_rf(i, mu_i, u_l) !! Get the reaction field contribution to the !! potential and torques: call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot) @@ -813,19 +892,19 @@ contains #endif - if (FF_uses_LJ .and. SimUsesLJ()) then + if (FF_uses_LJ .and. SIM_uses_LJ) then - if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) & + if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) & call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) endif - if (FF_uses_dipoles .and. SimUsesDipoles()) then + if (FF_uses_dipoles .and. SIM_uses_dipoles) then - if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then + if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, & do_pot, do_stress) - if (FF_uses_RF .and. SimUsesRF()) then + if (FF_uses_RF .and. SIM_uses_RF) then call accumulate_rf(i, j, r, u_l) call rf_correct_forces(i, j, d, r, u_l, f, do_stress) endif @@ -833,18 +912,18 @@ contains endif endif - if (FF_uses_Sticky .and. SimUsesSticky()) then + if (FF_uses_Sticky .and. SIM_uses_sticky) then - if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then + if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, & do_pot, do_stress) endif endif - if (FF_uses_GB .and. SimUsesGB()) then + if (FF_uses_GB .and. SIM_uses_GB) then - if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then + if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, & do_pot, do_stress) endif @@ -853,9 +932,9 @@ contains - if (FF_uses_EAM .and. SimUsesEAM()) then + if (FF_uses_EAM .and. SIM_uses_EAM) then - if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then + if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) endif @@ -900,14 +979,13 @@ contains #endif - if (FF_uses_EAM .and. SimUsesEAM()) then - call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i) - call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j) - - if ( is_EAM_i .and. is_EAM_j ) & + if (FF_uses_EAM .and. SIM_uses_EAM) then + + if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) & call calc_EAM_prepair_rho(i, j, d, r, rijsq ) - endif + endif + end subroutine do_prepair @@ -917,7 +995,7 @@ contains integer :: nlocal real( kind = dp ) :: pot - if (FF_uses_EAM .and. SimUsesEAM()) then + if (FF_uses_EAM .and. SIM_uses_EAM) then call calc_EAM_preforce_Frho(nlocal,pot) endif @@ -936,7 +1014,7 @@ contains d(1:3) = q_j(1:3) - q_i(1:3) ! Wrap back into periodic box if necessary - if ( SimUsesPBC() ) then + if ( SIM_uses_PBC ) then if( .not.boxIsOrthorhombic ) then ! calc the scaled coordinates. @@ -976,29 +1054,6 @@ contains end subroutine get_interatomic_vector - subroutine check_initialization(error) - integer, intent(out) :: error - - error = 0 - ! Make sure we are properly initialized. - if (.not. do_forces_initialized) then - write(*,*) "Forces not initialized" - error = -1 - return - endif - -#ifdef IS_MPI - if (.not. isMPISimSet()) then - write(default_error,*) "ERROR: mpiSimulation has not been initialized!" - error = -1 - return - endif -#endif - - return - end subroutine check_initialization - - subroutine zero_work_arrays() #ifdef IS_MPI @@ -1031,7 +1086,7 @@ contains #endif - if (FF_uses_EAM .and. SimUsesEAM()) then + if (FF_uses_EAM .and. SIM_uses_EAM) then call clean_EAM() endif