--- trunk/OOPSE/libmdtools/do_Forces.F90 2004/05/01 18:52:38 1144 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/05/26 16:41:23 1197 @@ -4,13 +4,14 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.53 2004-05-01 18:52:38 tim Exp $, $Date: 2004-05-01 18:52:38 $, $Name: not supported by cvs2svn $, $Revision: 1.53 $ +!! @version $Id: do_Forces.F90,v 1.62 2004-05-26 16:41:23 gezelter Exp $, $Date: 2004-05-26 16:41:23 $, $Name: not supported by cvs2svn $, $Revision: 1.62 $ module do_Forces use force_globals use simulation use definitions use atype_module + use switcheroo use neighborLists use lj use sticky_pair @@ -30,7 +31,11 @@ module do_Forces #define __FORTRAN90 #include "fForceField.h" +#include "fSwitchingFunction.h" + INTEGER, PARAMETER:: PREPAIR_LOOP = 1 + INTEGER, PARAMETER:: PAIR_LOOP = 2 + logical, save :: haveRlist = .false. logical, save :: haveNeighborList = .false. logical, save :: havePolicies = .false. @@ -63,7 +68,6 @@ module do_Forces public :: do_force_loop public :: setRlistDF - #ifdef PROFILE public :: getforcetime real, save :: forceTime = 0 @@ -279,13 +283,7 @@ contains !! Assume sanity (for the sake of argument) haveSaneForceField = .true. - !! - if (FF_uses_charges) then - dielect = getDielect() - call initialize_charge(dielect) - endif - !! check to make sure the FF_uses_RF setting makes sense if (FF_uses_dipoles) then @@ -364,20 +362,22 @@ contains endif haveNeighborList = .true. endif + + end subroutine init_FF !! Does force loop over i,j pairs. Calls do_pair to calculates forces. !-------------------------------------------------------------> - subroutine do_force_loop(q, qcom, A, u_l, f, t, tau, pot, & + subroutine do_force_loop(q, q_group, 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,nLocal) :: q + real ( kind = dp ), dimension(3, nLocal) :: q !! molecular center-of-mass position array - real ( kind = dp ), dimension(3,nLocal) :: qcom + real ( kind = dp ), dimension(3, nGroup) :: q_group !! Rotation Matrix for each long range particle in simulation. - real( kind = dp), dimension(9,nLocal) :: A + real( kind = dp), dimension(9, nLocal) :: A !! Unit vectors for dipoles (lab frame) real( kind = dp ), dimension(3,nLocal) :: u_l !! Force array provided by C, dimensioned by getNlocal @@ -391,25 +391,32 @@ contains logical ( kind = 2) :: do_pot_c, do_stress_c logical :: do_pot logical :: do_stress + logical :: in_switching_region #ifdef IS_MPI real( kind = DP ) :: pot_local integer :: nrow integer :: ncol integer :: nprocs + integer :: nrow_group + integer :: ncol_group #endif integer :: natoms logical :: update_nlist - integer :: i, j, jbeg, jend, jnab + integer :: i, j, jstart, jend, jnab + integer :: istart, iend + integer :: ia, jb, atom1, atom2 integer :: nlist - real( kind = DP ) :: rijsq, rcijsq - real(kind=dp),dimension(3) :: d, dc + real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij + real( kind = DP ) :: sw, dswdr, swderiv, mf + real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij real(kind=dp) :: rfpot, mu_i, virial - integer :: me_i, me_j + integer :: me_i, me_j, n_in_i, n_in_j logical :: is_dp_i integer :: neighborListSize integer :: listerror, error integer :: localError integer :: propPack_i, propPack_j + integer :: loopStart, loopEnd, loop real(kind=dp) :: listSkin = 1.0 @@ -419,6 +426,8 @@ contains pot_local = 0.0_dp nrow = getNrow(plan_row) ncol = getNcol(plan_col) + nrow_group = getNrowGroup(plan_row) + ncol_group = getNcolGroup(plan_col) #else natoms = nlocal #endif @@ -438,14 +447,12 @@ contains #ifdef IS_MPI - call gather(q,q_Row,plan_row3d) - call gather(q,q_Col,plan_col3d) - - if (SIM_uses_molecular_cutoffs) then - call gather(qcom,qcom_Row,plan_row3d) - call gather(qcom,qcom_Col,plan_col3d) - endif - + call gather(q, q_Row, plan_row3d) + call gather(q, q_Col, plan_col3d) + + call gather(q_group, q_group_Row, plan_row_Group_3d) + call gather(q_group, q_group_Col, plan_col_Group_3d) + 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) @@ -462,144 +469,74 @@ contains nloops = nloops + 1 #endif + loopEnd = PAIR_LOOP if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then - !! See if we need to update neighbor lists - if (SIM_uses_molecular_cutoffs) then - call checkNeighborList(nlocal, qcom, listSkin, update_nlist) - else - call checkNeighborList(nlocal, q, listSkin, update_nlist) + loopStart = PREPAIR_LOOP + else + loopStart = PAIR_LOOP + endif + + do loop = loopStart, loopEnd + + ! See if we need to update neighbor lists + ! (but only on the first time through): + if (loop .eq. loopStart) then + call checkNeighborList(nGroup, q_group, listSkin, update_nlist) endif - !! 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 - !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>> -#ifdef IS_MPI - if (update_nlist) then - - !! save current configuration, construct neighbor list, - !! and calculate forces - if (SIM_uses_molecular_cutoffs) then - call saveNeighborList(nlocal, qcom) - else - call saveNeighborList(nlocal, q) - endif - + !! save current configuration and construct neighbor list + call saveNeighborList(nGroup, q_group) 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 - - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), & - dc, rcijsq) - else - call get_interatomic_vector(q_Row(:,i), q_Col(:,j), & - dc, rcijsq) - endif - - if (rcijsq < 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 - - if (SIM_uses_molecular_cutoffs) then - ! since the simulation distances were in molecular COMs: - call get_interatomic_vector(q_Row(:,i), q_Col(:,j), & - d, rijsq) - else - d(1:3) = dc(1:3) - rijsq = rcijsq - endif - - call do_prepair(i, j, rijsq, d, rcijsq, dc, & - 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) - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom_Row(:,i),qcom_Col(:,j),& - dc, rcijsq) - else - dc(1:3) = d(1:3) - rcijsq = rijsq - endif - - call do_prepair(i, j, rijsq, d, rcijsq, dc, & - do_pot, do_stress, u_l, A, f, t, pot_local) - - enddo - endif - enddo + nlist = 0 endif + istart = 1 +#ifdef IS_MPI + iend = nrow_group #else - - if (update_nlist) then + iend = nGroup - 1 +#endif + outer: do i = istart, iend - ! save current configuration, contruct neighbor list, - ! and calculate forces - call saveNeighborList(natoms, q) + if (update_nlist) point(i) = nlist + 1 - neighborListSize = size(list) + n_in_i = groupStart(i+1) - groupStart(i) - nlist = 0 + if (update_nlist) then +#ifdef IS_MPI + jstart = 1 + jend = ncol_group +#else + jstart = i+1 + jend = nGroup +#endif + else + jstart = point(i) + jend = point(i+1) - 1 + ! make sure group i has neighbors + if (jstart .gt. jend) cycle outer + endif - do i = 1, natoms-1 - point(i) = nlist + 1 - - prepair_inner: do j = i+1, natoms - - if (skipThisPair(i,j)) cycle prepair_inner - - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom(:,i), qcom(:,j), & - dc, rcijsq) - else - call get_interatomic_vector(q(:,i), q(:,j), dc, rcijsq) - endif - - if (rcijsq < rlistsq) then - + do jnab = jstart, jend + if (update_nlist) then + j = jnab + else + j = list(jnab) + endif +#ifdef IS_MPI + call get_interatomic_vector(q_group_Row(:,i), & + q_group_Col(:,j), d_grp, rgrpsq) +#else + call get_interatomic_vector(q_group(:,i), & + q_group(:,j), d_grp, rgrpsq) +#endif + if (rgrpsq < rlistsq) then + if (update_nlist) then nlist = nlist + 1 if (nlist > neighborListSize) then - call expandNeighborList(natoms, listerror) + call expandNeighborList(nGroup, listerror) if (listerror /= 0) then error = -1 write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded." @@ -609,256 +546,134 @@ contains endif list(nlist) = j - - - if (SIM_uses_molecular_cutoffs) then - ! since the simulation distances were in molecular COMs: - call get_interatomic_vector(q(:,i), q(:,j), & - d, rijsq) - else - dc(1:3) = d(1:3) - rcijsq = rijsq - endif - - call do_prepair(i, j, rijsq, d, rcijsq, dc, & - 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) - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom(:,i), qcom(:,j), & - dc, rcijsq) - else - dc(1:3) = d(1:3) - rcijsq = rijsq - endif - - call do_prepair(i, j, rijsq, d, rcijsq, dc, & - 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) - - neighborListSize = size(list) - nlist = 0 - - do i = 1, nrow - - point(i) = nlist + 1 - - inner: do j = 1, ncol - - if (skipThisPair(i,j)) cycle inner - - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), & - dc, rcijsq) - else - call get_interatomic_vector(q_Row(:,i), q_Col(:,j), & - dc, rcijsq) - endif - - if (rcijsq < 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) + if (loop .eq. PAIR_LOOP) then + vij = 0.0d0 + fij(1:3) = 0.0d0 endif - list(nlist) = j + call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & + in_switching_region) - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(q_Row(:,i), q_Col(:,j), & - d, rijsq) - else - d(1:3) = dc(1:3) - rijsq = rcijsq - endif + n_in_j = groupStart(j+1) - groupStart(j) - call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, & - do_pot, do_stress, u_l, A, f, t, pot_local) - - endif - enddo 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) - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), & - dc, rcijsq) - else - dc(1:3) = d(1:3) - rcijsq = rijsq - endif - - call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, & - do_pot, do_stress, u_l, A, f, t, pot_local) - - enddo - endif - enddo - endif - + do ia = groupStart(i), groupStart(i+1)-1 + + atom1 = groupList(ia) + + inner: do jb = groupStart(j), groupStart(j+1)-1 + + atom2 = groupList(jb) + + if (skipThisPair(atom1, atom2)) cycle inner + + if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then + d_atm(1:3) = d_grp(1:3) + ratmsq = rgrpsq + else +#ifdef IS_MPI + call get_interatomic_vector(q_Row(:,atom1), & + q_Col(:,atom2), d_atm, ratmsq) #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 - - inner: do j = i+1, natoms - - if (skipThisPair(i,j)) cycle inner - - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom(:,i), qcom(:,j), & - dc, rcijsq) - else - call get_interatomic_vector(q(:,i), q(:,j), & - dc, rcijsq) - endif - - if (rcijsq < rlistsq) then + call get_interatomic_vector(q(:,atom1), & + q(:,atom2), d_atm, ratmsq) +#endif + endif + if (loop .eq. PREPAIR_LOOP) then +#ifdef IS_MPI + call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & + rgrpsq, d_grp, do_pot, do_stress, & + u_l, A, f, t, pot_local) +#else + call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & + rgrpsq, d_grp, do_pot, do_stress, & + u_l, A, f, t, pot) +#endif + else +#ifdef IS_MPI + call do_pair(atom1, atom2, ratmsq, d_atm, sw, & + do_pot, & + u_l, A, f, t, pot_local, vpair, fpair) +#else + call do_pair(atom1, atom2, ratmsq, d_atm, sw, & + do_pot, & + u_l, A, f, t, pot, vpair, fpair) +#endif + vij = vij + vpair + fij(1:3) = fij(1:3) + fpair(1:3) + endif + enddo inner + enddo - 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) + if (loop .eq. PAIR_LOOP) then + if (in_switching_region) then + swderiv = vij*dswdr/rgrp + fij(1) = fij(1) + swderiv*d_grp(1) + fij(2) = fij(2) + swderiv*d_grp(2) + fij(3) = fij(3) + swderiv*d_grp(3) + + do ia=groupStart(i), groupStart(i+1)-1 + atom1=groupList(ia) + mf = mfact(atom1) +#ifdef IS_MPI + f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf + f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf + f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf +#else + f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf + f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf + f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf +#endif + enddo + + do jb=groupStart(j), groupStart(j+1)-1 + atom2=groupList(jb) + mf = mfact(atom2) +#ifdef IS_MPI + f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf + f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf + f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf +#else + f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf + f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf + f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf +#endif + enddo + endif + + if (do_stress) call add_stress_tensor(d_grp, fij) endif - - list(nlist) = j - - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(q(:,i), q(:,j), & - d, rijsq) - else - d(1:3) = dc(1:3) - rijsq = rcijsq - endif - - call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, & - do_pot, do_stress, u_l, A, f, t, pot) - - endif - enddo inner - enddo + end if + enddo + enddo outer - 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) - if (SIM_uses_molecular_cutoffs) then - call get_interatomic_vector(qcom(:,i), qcom(:,j), & - dc, rcijsq) - else - dc(1:3) = d(1:3) - rcijsq = rijsq - endif - - call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, & - do_pot, do_stress, u_l, A, f, t, pot) - - enddo + if (update_nlist) then +#ifdef IS_MPI + point(nrow_group + 1) = nlist + 1 +#else + point(nGroup) = nlist + 1 +#endif + if (loop .eq. PREPAIR_LOOP) then + ! we just did the neighbor list update on the first + ! pass, so we don't need to do it + ! again on the second pass + update_nlist = .false. endif - enddo - endif + endif + + if (loop .eq. PREPAIR_LOOP) then + call do_preforce(nlocal, pot) + endif + + enddo -#endif - - ! phew, done with main loop. - !! Do timing #ifdef PROFILE call cpu_time(forceTimeFinal) forceTime = forceTime + forceTimeFinal - forceTimeInitial -#endif +#endif - #ifdef IS_MPI !!distribute forces @@ -904,7 +719,7 @@ contains do i = 1, nlocal pot_local = pot_local + pot_Temp(i) enddo - + endif #endif @@ -974,35 +789,31 @@ contains endif #endif - - + end subroutine do_force_loop - subroutine do_pair(i, j, rijsq, d, rcijsq, dc, mfact, do_pot, do_stress, & - u_l, A, f, t, pot) + subroutine do_pair(i, j, rijsq, d, sw, do_pot, & + u_l, A, f, t, pot, vpair, fpair) - real( kind = dp ) :: pot + real( kind = dp ) :: pot, vpair, sw + real( kind = dp ), dimension(3) :: fpair real( kind = dp ), dimension(nLocal) :: mfact 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 + logical, intent(inout) :: do_pot integer, intent(in) :: i, j - real ( kind = dp ), intent(inout) :: rijsq, rcijsq - real ( kind = dp ) :: r, rc - real ( kind = dp ), intent(inout) :: d(3), dc(3) + real ( kind = dp ), intent(inout) :: rijsq + real ( kind = dp ) :: r + real ( kind = dp ), intent(inout) :: d(3) integer :: me_i, me_j r = sqrt(rijsq) - if (SIM_uses_molecular_cutoffs) then - rc = sqrt(rcijsq) - else - rc = r - endif + vpair = 0.0d0 + fpair(1:3) = 0.0d0 - #ifdef IS_MPI if (tagRow(i) .eq. tagColumn(j)) then write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) @@ -1017,7 +828,13 @@ contains if (FF_uses_LJ .and. SIM_uses_LJ) 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) + !write(*,*) 'calling lj with' + !write(*,*) i, j, r, rijsq + !write(*,'(3es12.3)') d(1), d(2), d(3) + !write(*,'(3es12.3)') sw, vpair, pot + !write(*,*) + + call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) endif endif @@ -1025,8 +842,7 @@ contains 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, dc, rc, rcijsq, mfact, & - pot, f, do_pot, do_stress, SIM_uses_molecular_cutoffs) + call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) endif endif @@ -1034,11 +850,11 @@ contains 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) + call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, & + do_pot) 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) + call accumulate_rf(i, j, r, u_l, sw) + call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair) endif endif @@ -1047,8 +863,8 @@ contains if (FF_uses_Sticky .and. SIM_uses_sticky) 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) + call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, & + do_pot) endif endif @@ -1057,8 +873,8 @@ contains if (FF_uses_GB .and. SIM_uses_GB) 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) + call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, & + do_pot) endif endif @@ -1066,18 +882,18 @@ contains 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) + call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & + do_pot) endif endif - + end subroutine do_pair - - - subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, & + subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & do_pot, do_stress, u_l, A, f, t, pot) - real( kind = dp ) :: pot + + real( kind = dp ) :: pot, sw real( kind = dp ), dimension(3,nLocal) :: u_l real (kind=dp), dimension(9,nLocal) :: A real (kind=dp), dimension(3,nLocal) :: f @@ -1104,7 +920,7 @@ contains #ifdef IS_MPI if (tagRow(i) .eq. tagColumn(j)) then - write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) + write(0,*) 'do_prepair is doing', i , j, tagRow(i), tagColumn(j) endif me_i = atid_row(i) @@ -1116,19 +932,17 @@ contains 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 @@ -1199,8 +1013,8 @@ contains q_Row = 0.0_dp q_Col = 0.0_dp - qcom_Row = 0.0_dp - qcom_Col = 0.0_dp + q_group_Row = 0.0_dp + q_group_Col = 0.0_dp u_l_Row = 0.0_dp u_l_Col = 0.0_dp @@ -1307,19 +1121,10 @@ contains endif enddo - do i = 1, nExcludes_local - if (excludesLocal(1,i) == unique_id_1) then - if (excludesLocal(2,i) == unique_id_2) then - skip_it = .true. - return - endif - else - if (excludesLocal(1,i) == unique_id_2) then - if (excludesLocal(2,i) == unique_id_1) then - skip_it = .true. - return - endif - endif + do i = 1, nSkipsForAtom(unique_id_1) + if (skipsForAtom(unique_id_1, i) .eq. unique_id_2) then + skip_it = .true. + return endif end do @@ -1350,5 +1155,30 @@ contains #endif !! This cleans componets of force arrays belonging only to fortran + + subroutine add_stress_tensor(dpair, fpair) + + real( kind = dp ), dimension(3), intent(in) :: dpair, fpair + + ! because the d vector is the rj - ri vector, and + ! because fx, fy, fz are the force on atom i, we need a + ! negative sign here: + + tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1) + tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2) + tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3) + tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1) + tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2) + tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3) + tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1) + tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2) + tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3) + + !write(*,'(6es12.3)') fpair(1:3), tau_Temp(1), tau_Temp(5), tau_temp(9) + virial_Temp = virial_Temp + & + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) + + end subroutine add_stress_tensor end module do_Forces +