--- trunk/OOPSE/libmdtools/do_Forces.F90 2004/05/24 16:52:32 1191 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/05/24 21:03:30 1192 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.60 2004-05-21 15:58:40 gezelter Exp $, $Date: 2004-05-21 15:58:40 $, $Name: not supported by cvs2svn $, $Revision: 1.60 $ +!! @version $Id: do_Forces.F90,v 1.61 2004-05-24 21:03:25 gezelter Exp $, $Date: 2004-05-24 21:03:25 $, $Name: not supported by cvs2svn $, $Revision: 1.61 $ module do_Forces use force_globals @@ -400,11 +400,12 @@ contains integer :: natoms logical :: update_nlist integer :: i, j, jbeg, jend, jnab + integer :: istart, iend, jstart integer :: ia, jb, atom1, atom2 integer :: nlist real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij real( kind = DP ) :: sw, dswdr, swderiv, mf - real(kind=dp),dimension(3) :: d_atm, d_grp + real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij real(kind=dp) :: rfpot, mu_i, virial integer :: me_i, me_j, n_in_i, n_in_j logical :: is_dp_i @@ -475,26 +476,45 @@ contains !! 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 - + call saveNeighborList(nGroup, q_group) neighborListSize = size(list) nlist = 0 - do i = 1, nrow_group + istart = 1 +#ifdef IS_MPI + iend = nrow_group +#else + iend = nGroup - 1 +#endif + do i = istart, iend + point(i) = nlist + 1 - do j = 1, ncol_group + n_in_i = groupStart(i+1) - groupStart(i) + +#ifdef IS_MPI + jstart = 1 + jend = ncol_group +#else + jstart = i+1 + jend = nGroup +#endif + do j = jstart, jend +#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 nlist = nlist + 1 @@ -508,34 +528,70 @@ contains neighborListSize = size(list) endif - list(nlist) = j + list(nlist) = j + call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & + in_switching_region) + + n_in_j = groupStart(j+1) - groupStart(j) + do ia = groupStart(i), groupStart(i+1)-1 atom1 = groupList(ia) - + prepair_inner1: do jb = groupStart(j), groupStart(j+1)-1 - atom2 = groupList(jb) + atom2 = groupList(jb) if (skipThisPair(atom1, atom2)) cycle prepair_inner1 - call get_interatomic_vector(q_Row(:,atom1), & - q_Col(:,atom2), d_atm, ratmsq) - - call do_prepair(atom1, atom2, ratmsq, d_atm, & + 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 + call get_interatomic_vector(q(:,atom1), & + q(:,atom2), d_atm, ratmsq) +#endif + endif +#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 enddo prepair_inner1 enddo + end if enddo enddo - point(nrow_group + 1) = nlist + 1 - + +#ifdef IS_MPI + point(nrow_group + 1) = nlist + 1 +#else + point(nGroup) = nlist + 1 +#endif + else !! (of update_check) ! use the list to find the neighbors - do i = 1, nrow_group + + istart = 1 +#ifdef IS_MPI + iend = nrow_group +#else + iend = nGroup - 1 +#endif + + do i = istart, iend + + n_in_i = groupStart(i+1) - groupStart(i) + JBEG = POINT(i) JEND = POINT(i+1) - 1 ! check that group i has neighbors @@ -544,121 +600,55 @@ contains do jnab = jbeg, jend j = list(jnab) - do ia = groupStart(i), groupStart(i+1)-1 - atom1 = groupList(ia) - - prepair_inner2: do jb = groupStart(j), groupStart(j+1)-1 - atom2 = groupList(jb) - - if (skipThisPair(atom1, atom2)) cycle prepair_inner2 - - call get_interatomic_vector(q_Row(:,atom1), & - q_Col(:,atom2), d_atm, ratmsq) - - call do_prepair(atom1, atom2, ratmsq, d_atm, & - rgrpsq, d_grp, do_pot, do_stress, & - u_l, A, f, t, pot_local) - - enddo prepair_inner2 - enddo - enddo - endif - enddo - endif - +#ifdef IS_MPI + call get_interatomic_vector(q_group_Row(:,i), & + q_group_Col(:,j), d_grp, rgrpsq) #else - - if (update_nlist) then - - !! save current configuration, construct neighbor list, - !! and calculate forces - - call saveNeighborList(nGroup, q_group) - - neighborListSize = size(list) - nlist = 0 - - do i = 1, nGroup-1 - point(i) = nlist + 1 - - do j = i+1, nGroup - - call get_interatomic_vector(q_group(:,i), & - q_group(:,j), d_grp, rgrpsq) - - if (rgrpsq < rlistsq) then - nlist = nlist + 1 - - if (nlist > neighborListSize) then - call expandNeighborList(nGroup, 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 + call get_interatomic_vector(q_group(:,i), & + q_group(:,j), d_grp, rgrpsq) +#endif + call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & + in_switching_region) - list(nlist) = j + n_in_j = groupStart(j+1) - groupStart(j) do ia = groupStart(i), groupStart(i+1)-1 atom1 = groupList(ia) - prepair_inner1: do jb = groupStart(j), groupStart(j+1)-1 - atom2 = groupList(jb) - - if (skipThisPair(atom1, atom2)) cycle prepair_inner1 - - call get_interatomic_vector(q(:,atom1), & - q(:,atom2), d_atm, ratmsq) - - call do_prepair(atom1, atom2, ratmsq, d_atm, & - rgrpsq, d_grp, do_pot, do_stress, & - u_l, A, f, t, pot) - - enddo prepair_inner1 - enddo - end if - enddo - enddo - point(nGroup) = nlist + 1 - - else !! (of update_check) - - ! use the list to find the neighbors - do i = 1, nGroup-1 - JBEG = POINT(i) - JEND = POINT(i+1) - 1 - ! check that group i has neighbors - if (jbeg .le. jend) then - - do jnab = jbeg, jend - j = list(jnab) - - do ia = groupStart(i), groupStart(i+1)-1 - atom1 = groupList(ia) - prepair_inner2: do jb = groupStart(j), groupStart(j+1)-1 + atom2 = groupList(jb) if (skipThisPair(atom1, atom2)) cycle prepair_inner2 - - call get_interatomic_vector(q(:,atom1), & - q(:,atom2), d_atm, ratmsq) - call do_prepair(atom1, atom2, ratmsq, d_atm, & + 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 + call get_interatomic_vector(q(:,atom1), & + q(:,atom2), d_atm, ratmsq) +#endif + endif + +#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 enddo prepair_inner2 enddo enddo endif enddo endif - -#endif - !! Do rest of preforce calculations !! do necessary preforce calculations call do_preforce(nlocal,pot) @@ -667,12 +657,10 @@ contains else !! See if we need to update neighbor lists for non pre-pair call checkNeighborList(nGroup, q_group, listSkin, update_nlist) - endif + end if !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>> - -#ifdef IS_MPI - + if (update_nlist) then !! save current configuration, construct neighbor list, @@ -682,18 +670,35 @@ contains neighborListSize = size(list) nlist = 0 - - do i = 1, nrow_group + + istart = 1 +#ifdef IS_MPI + iend = nrow_group +#else + iend = nGroup - 1 +#endif + do i = istart, iend point(i) = nlist + 1 n_in_i = groupStart(i+1) - groupStart(i) - do j = 1, ncol_group +#ifdef IS_MPI + jstart = 1 + jend = ncol_group +#else + jstart = i+1 + jend = nGroup +#endif + do j = jstart, jend +#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 nlist = nlist + 1 @@ -709,8 +714,9 @@ contains list(nlist) = j - vij = 0.0d0 - + vij = 0.0d0 + fij(1:3) = 0.0d0 + call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & in_switching_region) @@ -728,49 +734,89 @@ contains 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 + call get_interatomic_vector(q(:,atom1), & + q(:,atom2), d_atm, ratmsq) +#endif endif - +#ifdef IS_MPI call do_pair(atom1, atom2, ratmsq, d_atm, sw, & - do_pot, do_stress, & - u_l, A, f, t, pot_local, vpair) - + 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) + enddo inner1 enddo 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) + end if enddo enddo - point(nrow_group + 1) = nlist + 1 +#ifdef IS_MPI + point(nrow_group + 1) = nlist + 1 +#else + point(nGroup) = nlist + 1 +#endif else !! (of update_check) ! use the list to find the neighbors - do i = 1, nrow_group + istart = 1 +#ifdef IS_MPI + iend = nrow_group +#else + iend = nGroup - 1 +#endif + + do i = istart, iend + n_in_i = groupStart(i+1) - groupStart(i) JBEG = POINT(i) @@ -781,10 +827,16 @@ contains do jnab = jbeg, jend j = list(jnab) +#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 vij = 0.0d0 + fij(1:3) = 0.0d0 + call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & in_switching_region) @@ -794,6 +846,7 @@ contains atom1 = groupList(ia) inner2: do jb = groupStart(j), groupStart(j+1)-1 + atom2 = groupList(jb) if (skipThisPair(atom1, atom2)) cycle inner2 @@ -802,213 +855,71 @@ contains 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 + call get_interatomic_vector(q(:,atom1), & + q(:,atom2), d_atm, ratmsq) +#endif endif - +#ifdef IS_MPI call do_pair(atom1, atom2, ratmsq, d_atm, sw, & - do_pot, do_stress, & - u_l, A, f, t, pot_local, vpair) - + 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) enddo inner2 enddo 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) + 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 - enddo - - do jb=groupStart(j), groupStart(j+1)-1 - atom2=groupList(jb) - mf = mfact(atom2) - 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 - enddo - endif - - enddo - endif - enddo - endif - #else - - if (update_nlist) then - - !! save current configuration, construct neighbor list, - !! and calculate forces - - call saveNeighborList(nGroup, q_group) - - neighborListSize = size(list) - nlist = 0 - - do i = 1, nGroup-1 - - point(i) = nlist + 1 - n_in_i = groupStart(i+1) - groupStart(i) - - do j = i+1, nGroup - - call get_interatomic_vector(q_group(:,i), & - q_group(:,j), d_grp, rgrpsq) - - if (rgrpsq < rlistsq) then - nlist = nlist + 1 - - if (nlist > neighborListSize) then - call expandNeighborList(nGroup, 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 - - vij = 0.0d0 - - call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & - in_switching_region) - - n_in_j = groupStart(j+1) - groupStart(j) - - do ia = groupStart(i), groupStart(i+1)-1 - atom1 = groupList(ia) - - inner1: do jb = groupStart(j), groupStart(j+1)-1 - atom2 = groupList(jb) - - if (skipThisPair(atom1, atom2)) cycle inner1 - - if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then - d_atm(1:3) = d_grp(1:3) - ratmsq = rgrpsq - else - call get_interatomic_vector(q(:,atom1), & - q(:,atom2), d_atm, ratmsq) - endif - - call do_pair(atom1, atom2, ratmsq, d_atm, sw, & - do_pot, do_stress, & - u_l, A, f, t, pot, vpair) - - vij = vij + vpair - - enddo inner1 - enddo - - if (in_switching_region) then - swderiv = vij*dswdr/rgrp - do ia=groupStart(i), groupStart(i+1)-1 - atom1=groupList(ia) - mf = mfact(atom1) 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 - end if - enddo - enddo - point(nGroup) = nlist + 1 - - else !! (of update_check) - - ! use the list to find the neighbors - do i = 1, nGroup-1 - - n_in_i = groupStart(i+1) - groupStart(i) - - JBEG = POINT(i) - JEND = POINT(i+1) - 1 - ! check that group i has neighbors - if (jbeg .le. jend) then - - do jnab = jbeg, jend - j = list(jnab) - - call get_interatomic_vector(q_group(:,i), & - q_group(:,j), d_grp, rgrpsq) + if (do_stress) call add_stress_tensor(d_grp, fij) - vij = 0.0d0 - - call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & - in_switching_region) - - n_in_j = groupStart(j+1) - groupStart(j) - - do ia = groupStart(i), groupStart(i+1)-1 - atom1 = groupList(ia) - - inner2: do jb = groupStart(j), groupStart(j+1)-1 - atom2 = groupList(jb) - - if (skipThisPair(atom1, atom2)) cycle inner2 - - if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then - d_atm(1:3) = d_grp(1:3) - ratmsq = rgrpsq - else - call get_interatomic_vector(q(:,atom1), & - q(:,atom2), d_atm, ratmsq) - endif - - call do_pair(atom1, atom2, ratmsq, d_atm, sw, & - do_pot, do_stress, & - u_l, A, f, t, pot, vpair) - - vij = vij + vpair - - enddo inner2 - enddo - - if (in_switching_region) then - swderiv = vij*dswdr/rgrp - - do ia=groupStart(i), groupStart(i+1)-1 - atom1=groupList(ia) - mf = mfact(atom1) - 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 - enddo - - do jb=groupStart(j), groupStart(j+1)-1 - atom2=groupList(jb) - mf = mfact(atom2) - 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 - enddo - endif enddo endif enddo endif -#endif - ! phew, done with main loop. !! Do timing @@ -1137,17 +1048,18 @@ contains end subroutine do_force_loop - subroutine do_pair(i, j, rijsq, d, sw, do_pot, do_stress, & - u_l, A, f, t, pot, vpair) + subroutine do_pair(i, j, rijsq, d, sw, do_pot, & + u_l, A, f, t, pot, vpair, fpair) 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 real ( kind = dp ) :: r @@ -1156,6 +1068,7 @@ contains r = sqrt(rijsq) vpair = 0.0d0 + fpair(1:3) = 0.0d0 #ifdef IS_MPI if (tagRow(i) .eq. tagColumn(j)) then @@ -1177,8 +1090,7 @@ contains !write(*,'(3es12.3)') sw, vpair, pot !write(*,*) - call do_lj_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, & - do_stress) + call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) endif endif @@ -1186,8 +1098,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, sw, vpair, pot, f, do_pot, & - do_stress) + call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) endif endif @@ -1195,11 +1106,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, sw, vpair, 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, sw) - call rf_correct_forces(i, j, d, r, u_l, sw, f, do_stress) + call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair) endif endif @@ -1208,8 +1119,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, sw, vpair, pot, A, 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 @@ -1218,8 +1129,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, sw, vpair, pot, u_l, 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 @@ -1227,17 +1138,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, sw, vpair, 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 @@ -1264,7 +1176,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) @@ -1276,19 +1188,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 @@ -1501,5 +1411,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 +