--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/07/01 22:39:53 571 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/04/15 16:18:26 1113 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.16 2003-07-01 22:39:53 gezelter Exp $, $Date: 2003-07-01 22:39:53 $, $Name: not supported by cvs2svn $, $Revision: 1.16 $ +!! @version $Id: do_Forces.F90,v 1.49 2004-04-15 16:18:26 tim Exp $, $Date: 2004-04-15 16:18:26 $, $Name: not supported by cvs2svn $, $Revision: 1.49 $ module do_Forces use force_globals @@ -15,8 +15,12 @@ module do_Forces use lj use sticky_pair use dipole_dipole + use charge_charge use reaction_field use gb_pair + use vector_class + use eam + use status #ifdef IS_MPI use mpiSimulation #endif @@ -27,19 +31,202 @@ module do_Forces #define __FORTRAN90 #include "fForceField.h" - logical, save :: do_forces_initialized = .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_charges 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_charges + 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 + public :: init_FF public :: do_force_loop + public :: setRlistDF + +#ifdef PROFILE + public :: getforcetime + real, save :: forceTime = 0 + real :: forceTimeInitial, forceTimeFinal + integer :: nLoops +#endif + + type :: Properties + logical :: is_lj = .false. + logical :: is_sticky = .false. + logical :: is_dp = .false. + logical :: is_gb = .false. + logical :: is_eam = .false. + logical :: is_charge = .false. + real(kind=DP) :: charge = 0.0_DP + real(kind=DP) :: dipole_moment = 0.0_DP + end type Properties + + type(Properties), dimension(:),allocatable :: PropertyMap + contains + subroutine setRlistDF( this_rlist ) + + real(kind=dp) :: this_rlist + + rlist = this_rlist + rlistsq = rlist * rlist + + haveRlist = .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_Charge", thisProperty) + PropertyMap(i)%is_Charge = thisProperty + + if (thisProperty) then + call getElementProperty(atypes, i, "charge", thisDPproperty) + PropertyMap(i)%charge = thisDPproperty + endif + + 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_charges = SimUsesCharges() + 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 @@ -64,13 +251,17 @@ contains FF_uses_LJ = .false. FF_uses_sticky = .false. + FF_uses_charges = .false. FF_uses_dipoles = .false. FF_uses_GB = .false. FF_uses_EAM = .false. call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList) if (nMatches .gt. 0) FF_uses_LJ = .true. - + + call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList) + if (nMatches .gt. 0) FF_uses_charges = .true. + call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList) if (nMatches .gt. 0) FF_uses_dipoles = .true. @@ -84,74 +275,88 @@ 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 - rrf = getRrf() - rt = getRt() - call initialize_dipole(rrf, rt) if (FF_uses_RF) then dielect = getDielect() - call initialize_rf(rrf, rt, dielect) + call initialize_rf(dielect) endif else if (FF_uses_RF) then write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' thisStat = -1 + haveSaneForceField = .false. return endif - endif + endif if (FF_uses_LJ) then - call getRcut(rcut) - select case (LJMIXPOLICY) case (LB_MIXING_RULE) - call init_lj_FF(LB_MIXING_RULE, rcut, my_status) + call init_lj_FF(LB_MIXING_RULE, my_status) case (EXPLICIT_MIXING_RULE) - call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status) + call init_lj_FF(EXPLICIT_MIXING_RULE, my_status) 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 - + + + if (FF_uses_EAM) then + call init_EAM_FF(my_status) + if (my_status /= 0) then + 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(getNlocal(), my_status) + call expandNeighborList(nLocal, my_status) if (my_Status /= 0) then write(default_error,*) "SimSetup: ExpandNeighborList returned error." thisStat = -1 return endif + haveNeighborList = .true. endif - - do_forces_initialized = .true. - + end subroutine init_FF @@ -160,15 +365,16 @@ contains subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, & error) !! Position array provided by C, dimensioned by getNlocal - real ( kind = dp ), dimension(3,getNlocal()) :: q + real ( kind = dp ), dimension(3,nLocal) :: q !! Rotation Matrix for each long range particle in simulation. - real( kind = dp), dimension(9,getNlocal()) :: A + real( kind = dp), dimension(9,nLocal) :: A !! Unit vectors for dipoles (lab frame) - real( kind = dp ), dimension(3,getNlocal()) :: u_l + real( kind = dp ), dimension(3,nLocal) :: u_l !! Force array provided by C, dimensioned by getNlocal - real ( kind = dp ), dimension(3,getNlocal()) :: f + real ( kind = dp ), dimension(3,nLocal) :: f !! Torsion array provided by C, dimensioned by getNlocal - real( kind = dp ), dimension(3,getNlocal()) :: t + real( kind = dp ), dimension(3,nLocal) :: t + !! Stress Tensor real( kind = dp), dimension(9) :: tau real ( kind = dp ) :: pot @@ -179,38 +385,37 @@ contains real( kind = DP ) :: pot_local integer :: nrow integer :: ncol + integer :: nprocs #endif - integer :: nlocal integer :: natoms logical :: update_nlist integer :: i, j, jbeg, jend, jnab integer :: nlist - real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut + real( kind = DP ) :: rijsq real(kind=dp),dimension(3) :: d real(kind=dp) :: rfpot, mu_i, virial - integer :: me_i + integer :: me_i, me_j logical :: is_dp_i integer :: neighborListSize integer :: listerror, error integer :: localError + integer :: propPack_i, propPack_j + real(kind=dp) :: listSkin = 1.0 + !! initialize local variables #ifdef IS_MPI pot_local = 0.0_dp - nlocal = getNlocal() nrow = getNrow(plan_row) ncol = getNcol(plan_col) #else - nlocal = getNlocal() natoms = nlocal #endif - - call getRcut(rcut,rc2=rcutsq) - call getRlist(rlist,rlistsq) - - call check_initialization(localError) + + call doReadyCheck(localError) if ( localError .ne. 0 ) then + call handleError("do_force_loop", "Not Initialized") error = -1 return end if @@ -226,7 +431,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) @@ -235,23 +440,177 @@ contains endif #endif - - if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then + +!! Begin force loop timing: +#ifdef PROFILE + call cpu_time(forceTimeInitial) + nloops = nloops + 1 +#endif + + if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then !! See if we need to update neighbor lists - call checkNeighborList(nlocal, q, rcut, rlist, update_nlist) + call checkNeighborList(nlocal, q, listSkin, update_nlist) !! if_mpi_gather_stuff_for_prepair !! do_prepair_loop_if_needed !! if_mpi_scatter_stuff_from_prepair !! if_mpi_gather_stuff_from_prepair_to_main_loop - else - !! See if we need to update neighbor lists - call checkNeighborList(nlocal, q, rcut, rlist, update_nlist) + +!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>> +#ifdef IS_MPI + + if (update_nlist) then + + !! save current configuration, construct neighbor list, + !! and calculate forces + call saveNeighborList(nlocal, q) + + neighborListSize = size(list) + nlist = 0 + + do i = 1, nrow + point(i) = nlist + 1 + + prepair_inner: do j = 1, ncol + + if (skipThisPair(i,j)) cycle prepair_inner + + call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) + + if (rijsq < rlistsq) then + + nlist = nlist + 1 + + if (nlist > neighborListSize) then + call expandNeighborList(nlocal, listerror) + if (listerror /= 0) then + error = -1 + write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded." + return + end if + neighborListSize = size(list) + endif + + list(nlist) = j + call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local) + endif + enddo prepair_inner + enddo + + point(nrow + 1) = nlist + 1 + + else !! (of update_check) + + ! use the list to find the neighbors + do i = 1, nrow + JBEG = POINT(i) + JEND = POINT(i+1) - 1 + ! check thiat molecule i has neighbors + if (jbeg .le. jend) then + + do jnab = jbeg, jend + j = list(jnab) + + call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) + call do_prepair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot_local) + + enddo + endif + enddo endif -#ifdef IS_MPI +#else if (update_nlist) then + ! save current configuration, contruct neighbor list, + ! and calculate forces + call saveNeighborList(natoms, q) + + neighborListSize = size(list) + + nlist = 0 + + do i = 1, natoms-1 + point(i) = nlist + 1 + + prepair_inner: do j = i+1, natoms + + if (skipThisPair(i,j)) cycle prepair_inner + + call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) + + + if (rijsq < rlistsq) then + + + nlist = nlist + 1 + + if (nlist > neighborListSize) then + call expandNeighborList(natoms, listerror) + if (listerror /= 0) then + error = -1 + write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded." + return + end if + neighborListSize = size(list) + endif + + list(nlist) = j + + call do_prepair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot) + + endif + enddo prepair_inner + enddo + + point(natoms) = nlist + 1 + + else !! (update) + + ! use the list to find the neighbors + do i = 1, natoms-1 + JBEG = POINT(i) + JEND = POINT(i+1) - 1 + ! check thiat molecule i has neighbors + if (jbeg .le. jend) then + + do jnab = jbeg, jend + j = list(jnab) + + call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) + call do_prepair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot) + + enddo + endif + enddo + endif +#endif + !! Do rest of preforce calculations + !! do necessary preforce calculations + call do_preforce(nlocal,pot) + ! we have already updated the neighbor list set it to false... + update_nlist = .false. + else + !! See if we need to update neighbor lists for non pre-pair + call checkNeighborList(nlocal, q, listSkin, update_nlist) + endif + + + + + +!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>> + + + + + +#ifdef IS_MPI + + if (update_nlist) then !! save current configuration, construct neighbor list, !! and calculate forces call saveNeighborList(nlocal, q) @@ -260,6 +619,7 @@ contains nlist = 0 do i = 1, nrow + point(i) = nlist + 1 inner: do j = 1, ncol @@ -268,7 +628,7 @@ contains call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) - if (rijsq < rlistsq) then + if (rijsq < rlistsq) then nlist = nlist + 1 @@ -284,10 +644,9 @@ contains list(nlist) = j - if (rijsq < rcutsq) then - call do_pair(i, j, rijsq, d, do_pot, do_stress, & - u_l, A, f, t, pot_local) - endif + call do_pair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot_local) + endif enddo inner enddo @@ -318,7 +677,7 @@ contains #else if (update_nlist) then - + ! save current configuration, contruct neighbor list, ! and calculate forces call saveNeighborList(natoms, q) @@ -337,7 +696,7 @@ contains call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) - if (rijsq < rlistsq) then + if (rijsq < rlistsq) then nlist = nlist + 1 @@ -353,10 +712,9 @@ contains list(nlist) = j - if (rijsq < rcutsq) then - call do_pair(i, j, rijsq, d, do_pot, do_stress, & + call do_pair(i, j, rijsq, d, do_pot, do_stress, & u_l, A, f, t, pot) - endif + endif enddo inner enddo @@ -387,7 +745,14 @@ contains #endif ! phew, done with main loop. - + +!! Do timing +#ifdef PROFILE + call cpu_time(forceTimeFinal) + forceTime = forceTime + forceTimeFinal - forceTimeInitial +#endif + + #ifdef IS_MPI !!distribute forces @@ -403,7 +768,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 @@ -437,10 +802,10 @@ 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) call scatter(rf_Col,rf_Temp,plan_col3d) @@ -449,7 +814,7 @@ contains end do #endif - do i = 1, getNlocal() + do i = 1, nLocal rfpot = 0.0_DP #ifdef IS_MPI @@ -457,12 +822,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) @@ -499,99 +866,160 @@ contains tau = tau_Temp virial = virial_Temp endif - + #endif + + end subroutine do_force_loop subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) real( kind = dp ) :: pot - real( kind = dp ), dimension(3,getNlocal()) :: u_l - real (kind=dp), dimension(9,getNlocal()) :: A - real (kind=dp), dimension(3,getNlocal()) :: f - real (kind=dp), dimension(3,getNlocal()) :: t + real( kind = dp ), dimension(3,nLocal) :: u_l + real (kind=dp), dimension(9,nLocal) :: A + real (kind=dp), dimension(3,nLocal) :: f + real (kind=dp), dimension(3,nLocal) :: t logical, intent(inout) :: do_pot, do_stress integer, intent(in) :: i, j real ( kind = dp ), intent(inout) :: rijsq real ( kind = dp ) :: r real ( kind = dp ), intent(inout) :: d(3) - logical :: is_LJ_i, is_LJ_j - logical :: is_DP_i, is_DP_j - logical :: is_GB_i, is_GB_j - logical :: is_Sticky_i, is_Sticky_j integer :: me_i, me_j r = sqrt(rijsq) - - #ifdef IS_MPI if (tagRow(i) .eq. tagColumn(j)) then write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) endif - me_i = atid_row(i) me_j = atid_col(j) - #else - me_i = atid(i) me_j = atid(j) - #endif - - if (FF_uses_LJ .and. SimUsesLJ()) then - call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i) - call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j) - - if ( is_LJ_i .and. is_LJ_j ) & - call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) - endif - - if (FF_uses_dipoles .and. SimUsesDipoles()) then - call getElementProperty(atypes, me_i, "is_DP", is_DP_i) - call getElementProperty(atypes, me_j, "is_DP", is_DP_j) + + if (FF_uses_LJ .and. SIM_uses_LJ) then - if ( is_DP_i .and. is_DP_j ) then - + if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then + call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) + endif + + endif + + if (FF_uses_charges .and. SIM_uses_charges) then + + if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then + call do_charge_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) + endif + + endif + + if (FF_uses_dipoles .and. SIM_uses_dipoles) 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 - + endif endif + endif - if (FF_uses_Sticky .and. SimUsesSticky()) then + if (FF_uses_Sticky .and. SIM_uses_sticky) then - call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i) - call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j) - - if ( is_Sticky_i .and. is_Sticky_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 - - call getElementProperty(atypes, me_i, "is_GB", is_GB_i) - call getElementProperty(atypes, me_j, "is_GB", is_GB_j) + if (FF_uses_GB .and. SIM_uses_GB) then - if ( is_GB_i .and. is_GB_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 + endif - + + if (FF_uses_EAM .and. SIM_uses_EAM) 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 + + endif + end subroutine do_pair + + + + subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) + real( kind = dp ) :: pot + real( kind = dp ), dimension(3,nLocal) :: u_l + real (kind=dp), dimension(9,nLocal) :: A + real (kind=dp), dimension(3,nLocal) :: f + real (kind=dp), dimension(3,nLocal) :: t + + logical, intent(inout) :: do_pot, do_stress + integer, intent(in) :: i, j + real ( kind = dp ), intent(inout) :: rijsq + real ( kind = dp ) :: r + real ( kind = dp ), intent(inout) :: d(3) + + logical :: is_EAM_i, is_EAM_j + + integer :: me_i, me_j + + r = sqrt(rijsq) + + +#ifdef IS_MPI + if (tagRow(i) .eq. tagColumn(j)) then + write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) + endif + + me_i = atid_row(i) + me_j = atid_col(j) + +#else + + me_i = atid(i) + me_j = atid(j) + +#endif + + 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 + + end subroutine do_prepair + + + + + subroutine do_preforce(nlocal,pot) + integer :: nlocal + real( kind = dp ) :: pot + + if (FF_uses_EAM .and. SIM_uses_EAM) then + call calc_EAM_preforce_Frho(nlocal,pot) + endif + + + end subroutine do_preforce + + subroutine get_interatomic_vector(q_i, q_j, d, r_sq) real (kind = dp), dimension(3) :: q_i @@ -603,23 +1031,22 @@ 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. - scaled = matmul(d, HmatInv) + scaled = matmul(HmatInv, d) ! wrap the scaled coordinates - do i = 1, 3 - scaled(i) = scaled(i) - anint(scaled(i)) - enddo + scaled = scaled - anint(scaled) + ! calc the wrapped real coordinates from the wrapped scaled ! coordinates - d = matmul(scaled,Hmat) + d = matmul(Hmat,scaled) else ! calc the scaled coordinates. @@ -644,28 +1071,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 - 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 @@ -697,6 +1102,15 @@ contains #endif + + if (FF_uses_EAM .and. SIM_uses_EAM) then + call clean_EAM() + endif + + + + + rf = 0.0_dp tau_Temp = 0.0_dp virial_Temp = 0.0_dp @@ -809,4 +1223,13 @@ end module do_Forces doesit = FF_uses_RF end function FF_RequiresPostpairCalc +#ifdef PROFILE + function getforcetime() result(totalforcetime) + real(kind=dp) :: totalforcetime + totalforcetime = forcetime + end function getforcetime +#endif + +!! This cleans componets of force arrays belonging only to fortran + end module do_Forces