45 |
|
|
46 |
|
!! @author Charles F. Vardeman II |
47 |
|
!! @author Matthew Meineke |
48 |
< |
!! @version $Id: doForces.F90,v 1.85 2007-04-06 21:53:41 gezelter Exp $, $Date: 2007-04-06 21:53:41 $, $Name: not supported by cvs2svn $, $Revision: 1.85 $ |
48 |
> |
!! @version $Id: doForces.F90,v 1.93 2008-05-27 16:39:05 chuckv Exp $, $Date: 2008-05-27 16:39:05 $, $Name: not supported by cvs2svn $, $Revision: 1.93 $ |
49 |
|
|
50 |
|
|
51 |
|
module doForces |
62 |
|
use shapes |
63 |
|
use vector_class |
64 |
|
use eam |
65 |
+ |
use MetalNonMetal |
66 |
|
use suttonchen |
67 |
|
use status |
68 |
|
#ifdef IS_MPI |
97 |
|
logical, save :: FF_uses_GayBerne |
98 |
|
logical, save :: FF_uses_EAM |
99 |
|
logical, save :: FF_uses_SC |
100 |
< |
logical, save :: FF_uses_MEAM |
100 |
> |
logical, save :: FF_uses_MNM |
101 |
|
|
102 |
|
|
103 |
|
logical, save :: SIM_uses_DirectionalAtoms |
104 |
|
logical, save :: SIM_uses_EAM |
105 |
|
logical, save :: SIM_uses_SC |
106 |
< |
logical, save :: SIM_uses_MEAM |
106 |
> |
logical, save :: SIM_uses_MNM |
107 |
|
logical, save :: SIM_requires_postpair_calc |
108 |
|
logical, save :: SIM_requires_prepair_calc |
109 |
|
logical, save :: SIM_uses_PBC |
114 |
|
|
115 |
|
real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut |
116 |
|
real(kind=dp), save :: skinThickness |
117 |
< |
logical, save :: defaultDoShift |
117 |
> |
logical, save :: defaultDoShiftPot |
118 |
> |
logical, save :: defaultDoShiftFrc |
119 |
|
|
120 |
|
public :: init_FF |
121 |
|
public :: setCutoffs |
171 |
|
logical :: i_is_EAM |
172 |
|
logical :: i_is_Shape |
173 |
|
logical :: i_is_SC |
172 |
– |
logical :: i_is_MEAM |
174 |
|
logical :: j_is_LJ |
175 |
|
logical :: j_is_Elect |
176 |
|
logical :: j_is_Sticky |
179 |
|
logical :: j_is_EAM |
180 |
|
logical :: j_is_Shape |
181 |
|
logical :: j_is_SC |
181 |
– |
logical :: j_is_MEAM |
182 |
|
real(kind=dp) :: myRcut |
183 |
|
|
184 |
|
if (.not. associated(atypes)) then |
216 |
|
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
217 |
|
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
218 |
|
call getElementProperty(atypes, i, "is_SC", i_is_SC) |
219 |
– |
call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM) |
219 |
|
|
220 |
|
do j = i, nAtypes |
221 |
|
|
230 |
|
call getElementProperty(atypes, j, "is_EAM", j_is_EAM) |
231 |
|
call getElementProperty(atypes, j, "is_Shape", j_is_Shape) |
232 |
|
call getElementProperty(atypes, j, "is_SC", j_is_SC) |
234 |
– |
call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM) |
233 |
|
|
234 |
|
if (i_is_LJ .and. j_is_LJ) then |
235 |
|
iHash = ior(iHash, LJ_PAIR) |
258 |
|
if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR) |
259 |
|
if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ) |
260 |
|
if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ) |
261 |
+ |
|
262 |
+ |
if ((i_is_EAM.or.i_is_SC).and.(.not.(j_is_EAM.or.j_is_SC))) iHash = ior(iHash, MNM_PAIR) |
263 |
+ |
if ((j_is_EAM.or.j_is_SC).and.(.not.(i_is_EAM.or.i_is_SC))) iHash = ior(iHash, MNM_PAIR) |
264 |
|
|
265 |
|
if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR) |
266 |
|
if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ) |
563 |
|
haveGtypeCutoffMap = .true. |
564 |
|
end subroutine createGtypeCutoffMap |
565 |
|
|
566 |
< |
subroutine setCutoffs(defRcut, defRsw) |
566 |
> |
subroutine setCutoffs(defRcut, defRsw, defSP, defSF) |
567 |
|
|
568 |
|
real(kind=dp),intent(in) :: defRcut, defRsw |
569 |
+ |
logical, intent(in) :: defSP, defSF |
570 |
|
character(len = statusMsgSize) :: errMsg |
571 |
|
integer :: localError |
572 |
|
|
573 |
|
defaultRcut = defRcut |
574 |
|
defaultRsw = defRsw |
575 |
|
|
576 |
< |
defaultDoShift = .false. |
576 |
> |
defaultDoShiftPot = defSP |
577 |
> |
defaultDoShiftFrc = defSF |
578 |
> |
|
579 |
|
if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then |
580 |
< |
|
581 |
< |
write(errMsg, *) & |
582 |
< |
'cutoffRadius and switchingRadius are set to the same', newline & |
583 |
< |
// tab, 'value. OOPSE will use shifted ', newline & |
584 |
< |
// tab, 'potentials instead of switching functions.' |
585 |
< |
|
586 |
< |
call handleInfo("setCutoffs", errMsg) |
587 |
< |
|
588 |
< |
defaultDoShift = .true. |
589 |
< |
|
580 |
> |
if (defaultDoShiftFrc) then |
581 |
> |
write(errMsg, *) & |
582 |
> |
'cutoffRadius and switchingRadius are set to the', newline & |
583 |
> |
// tab, 'same value. OOPSE will use shifted force', newline & |
584 |
> |
// tab, 'potentials instead of switching functions.' |
585 |
> |
|
586 |
> |
call handleInfo("setCutoffs", errMsg) |
587 |
> |
else |
588 |
> |
write(errMsg, *) & |
589 |
> |
'cutoffRadius and switchingRadius are set to the', newline & |
590 |
> |
// tab, 'same value. OOPSE will use shifted', newline & |
591 |
> |
// tab, 'potentials instead of switching functions.' |
592 |
> |
|
593 |
> |
call handleInfo("setCutoffs", errMsg) |
594 |
> |
|
595 |
> |
defaultDoShiftPot = .true. |
596 |
> |
endif |
597 |
> |
|
598 |
|
endif |
599 |
|
|
600 |
|
localError = 0 |
601 |
< |
call setLJDefaultCutoff( defaultRcut, defaultDoShift ) |
601 |
> |
call setLJDefaultCutoff( defaultRcut, defaultDoShiftPot, & |
602 |
> |
defaultDoShiftFrc ) |
603 |
|
call setElectrostaticCutoffRadius( defaultRcut, defaultRsw ) |
604 |
|
call setCutoffEAM( defaultRcut ) |
605 |
|
call setCutoffSC( defaultRcut ) |
606 |
+ |
call setMnMDefaultCutoff( defaultRcut, defaultDoShiftPot, & |
607 |
+ |
defaultDoShiftFrc ) |
608 |
|
call set_switch(defaultRsw, defaultRcut) |
609 |
|
call setHmatDangerousRcutValue(defaultRcut) |
610 |
|
|
794 |
|
|
795 |
|
!! Does force loop over i,j pairs. Calls do_pair to calculates forces. |
796 |
|
!-------------------------------------------------------------> |
797 |
< |
subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, & |
797 |
> |
subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, particle_pot,& |
798 |
|
do_pot_c, do_stress_c, error) |
799 |
|
!! Position array provided by C, dimensioned by getNlocal |
800 |
|
real ( kind = dp ), dimension(3, nLocal) :: q |
812 |
|
!! Stress Tensor |
813 |
|
real( kind = dp), dimension(9) :: tau |
814 |
|
real ( kind = dp ),dimension(LR_POT_TYPES) :: pot |
815 |
+ |
real( kind = dp ), dimension(nLocal) :: particle_pot |
816 |
|
logical ( kind = 2) :: do_pot_c, do_stress_c |
817 |
|
logical :: do_pot |
818 |
|
logical :: do_stress |
1095 |
|
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
1096 |
|
do_pot, eFrame, A, f, t, pot_local, vpair, & |
1097 |
|
fpair, d_grp, rgrp, rCut) |
1098 |
+ |
! particle_pot will be accumulated from row & column |
1099 |
+ |
! arrays later |
1100 |
|
#else |
1101 |
|
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
1102 |
< |
do_pot, eFrame, A, f, t, pot, vpair, fpair, & |
1103 |
< |
d_grp, rgrp, rCut) |
1102 |
> |
do_pot, eFrame, A, f, t, pot, vpair, & |
1103 |
> |
fpair, d_grp, rgrp, rCut) |
1104 |
> |
particle_pot(atom1) = particle_pot(atom1) + vpair*sw |
1105 |
> |
particle_pot(atom2) = particle_pot(atom2) + vpair*sw |
1106 |
|
#endif |
1107 |
|
vij = vij + vpair |
1108 |
|
fij(1) = fij(1) + fpair(1) |
1109 |
|
fij(2) = fij(2) + fpair(2) |
1110 |
|
fij(3) = fij(3) + fpair(3) |
1111 |
< |
if (do_stress.and.SIM_uses_AtomicVirial) then |
1111 |
> |
if (do_stress) then |
1112 |
|
call add_stress_tensor(d_atm, fpair, tau) |
1113 |
|
endif |
1114 |
|
endif |
1118 |
|
if (loop .eq. PAIR_LOOP) then |
1119 |
|
if (in_switching_region) then |
1120 |
|
swderiv = vij*dswdr/rgrp |
1121 |
< |
fij(1) = fij(1) + swderiv*d_grp(1) |
1122 |
< |
fij(2) = fij(2) + swderiv*d_grp(2) |
1123 |
< |
fij(3) = fij(3) + swderiv*d_grp(3) |
1121 |
> |
fg = swderiv*d_grp |
1122 |
> |
|
1123 |
> |
fij(1) = fij(1) + fg(1) |
1124 |
> |
fij(2) = fij(2) + fg(2) |
1125 |
> |
fij(3) = fij(3) + fg(3) |
1126 |
|
|
1127 |
+ |
if (do_stress .and. (n_in_i .eq. 1).and.(n_in_j .eq. 1)) then |
1128 |
+ |
call add_stress_tensor(d_atm, fg, tau) |
1129 |
+ |
endif |
1130 |
+ |
|
1131 |
|
do ia=groupStartRow(i), groupStartRow(i+1)-1 |
1132 |
|
atom1=groupListRow(ia) |
1133 |
|
mf = mfactRow(atom1) |
1143 |
|
f(2,atom1) = f(2,atom1) + fg(2) |
1144 |
|
f(3,atom1) = f(3,atom1) + fg(3) |
1145 |
|
#endif |
1146 |
< |
if (do_stress.and.SIM_uses_AtomicVirial) then |
1147 |
< |
! find the distance between the atom and the center of |
1148 |
< |
! the cutoff group: |
1146 |
> |
if (n_in_i .gt. 1) then |
1147 |
> |
if (do_stress.and.SIM_uses_AtomicVirial) then |
1148 |
> |
! find the distance between the atom and the center of |
1149 |
> |
! the cutoff group: |
1150 |
|
#ifdef IS_MPI |
1151 |
< |
call get_interatomic_vector(q_Row(:,atom1), & |
1152 |
< |
q_group_Row(:,i), dag, rag) |
1151 |
> |
call get_interatomic_vector(q_Row(:,atom1), & |
1152 |
> |
q_group_Row(:,i), dag, rag) |
1153 |
|
#else |
1154 |
< |
call get_interatomic_vector(q(:,atom1), & |
1155 |
< |
q_group(:,i), dag, rag) |
1154 |
> |
call get_interatomic_vector(q(:,atom1), & |
1155 |
> |
q_group(:,i), dag, rag) |
1156 |
|
#endif |
1157 |
< |
call add_stress_tensor(dag,fg,tau) |
1157 |
> |
call add_stress_tensor(dag,fg,tau) |
1158 |
> |
endif |
1159 |
|
endif |
1160 |
|
enddo |
1161 |
|
|
1174 |
|
f(2,atom2) = f(2,atom2) + fg(2) |
1175 |
|
f(3,atom2) = f(3,atom2) + fg(3) |
1176 |
|
#endif |
1177 |
< |
if (do_stress.and.SIM_uses_AtomicVirial) then |
1178 |
< |
! find the distance between the atom and the center of |
1179 |
< |
! the cutoff group: |
1177 |
> |
if (n_in_j .gt. 1) then |
1178 |
> |
if (do_stress.and.SIM_uses_AtomicVirial) then |
1179 |
> |
! find the distance between the atom and the center of |
1180 |
> |
! the cutoff group: |
1181 |
|
#ifdef IS_MPI |
1182 |
< |
call get_interatomic_vector(q_Col(:,atom2), & |
1183 |
< |
q_group_Col(:,j), dag, rag) |
1182 |
> |
call get_interatomic_vector(q_Col(:,atom2), & |
1183 |
> |
q_group_Col(:,j), dag, rag) |
1184 |
|
#else |
1185 |
< |
call get_interatomic_vector(q(:,atom2), & |
1186 |
< |
q_group(:,j), dag, rag) |
1185 |
> |
call get_interatomic_vector(q(:,atom2), & |
1186 |
> |
q_group(:,j), dag, rag) |
1187 |
|
#endif |
1188 |
< |
call add_stress_tensor(dag,fg,tau) |
1189 |
< |
endif |
1190 |
< |
|
1188 |
> |
call add_stress_tensor(dag,fg,tau) |
1189 |
> |
endif |
1190 |
> |
endif |
1191 |
|
enddo |
1192 |
|
endif |
1193 |
< |
if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then |
1194 |
< |
call add_stress_tensor(d_grp, fij, tau) |
1195 |
< |
endif |
1193 |
> |
!if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then |
1194 |
> |
! call add_stress_tensor(d_grp, fij, tau) |
1195 |
> |
!endif |
1196 |
|
endif |
1197 |
|
endif |
1198 |
|
endif |
1215 |
|
endif |
1216 |
|
|
1217 |
|
if (loop .eq. PREPAIR_LOOP) then |
1218 |
+ |
#ifdef IS_MPI |
1219 |
+ |
call do_preforce(nlocal, pot_local) |
1220 |
+ |
#else |
1221 |
|
call do_preforce(nlocal, pot) |
1222 |
+ |
#endif |
1223 |
|
endif |
1224 |
|
|
1225 |
|
enddo |
1271 |
|
+ pot_Temp(1:LR_POT_TYPES,i) |
1272 |
|
enddo |
1273 |
|
|
1274 |
+ |
do i = 1,LR_POT_TYPES |
1275 |
+ |
particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal) |
1276 |
+ |
enddo |
1277 |
+ |
|
1278 |
|
pot_Temp = 0.0_DP |
1279 |
+ |
|
1280 |
|
do i = 1,LR_POT_TYPES |
1281 |
|
call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col) |
1282 |
|
end do |
1283 |
+ |
|
1284 |
|
do i = 1, nlocal |
1285 |
|
pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)& |
1286 |
|
+ pot_Temp(1:LR_POT_TYPES,i) |
1287 |
|
enddo |
1288 |
|
|
1289 |
+ |
do i = 1,LR_POT_TYPES |
1290 |
+ |
particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal) |
1291 |
+ |
enddo |
1292 |
+ |
|
1293 |
+ |
|
1294 |
|
endif |
1295 |
|
#endif |
1296 |
|
|
1441 |
|
|
1442 |
|
real( kind = dp ) :: vpair, sw |
1443 |
|
real( kind = dp ), dimension(LR_POT_TYPES) :: pot |
1444 |
+ |
real( kind = dp ), dimension(nLocal) :: particle_pot |
1445 |
|
real( kind = dp ), dimension(3) :: fpair |
1446 |
|
real( kind = dp ), dimension(nLocal) :: mfact |
1447 |
|
real( kind = dp ), dimension(9,nLocal) :: eFrame |
1456 |
|
real ( kind = dp ), intent(inout) :: d(3) |
1457 |
|
real ( kind = dp ), intent(inout) :: d_grp(3) |
1458 |
|
real ( kind = dp ), intent(inout) :: rCut |
1459 |
< |
real ( kind = dp ) :: r |
1459 |
> |
real ( kind = dp ) :: r, pair_pot |
1460 |
|
real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx |
1461 |
|
integer :: me_i, me_j |
1462 |
|
integer :: k |
1528 |
|
pot(METALLIC_POT), f, do_pot) |
1529 |
|
endif |
1530 |
|
|
1531 |
+ |
if ( iand(iHash, MNM_PAIR).ne.0 ) then |
1532 |
+ |
call do_mnm_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & |
1533 |
+ |
pot(VDW_POT), A, f, t, do_pot) |
1534 |
+ |
endif |
1535 |
+ |
|
1536 |
|
end subroutine do_pair |
1537 |
|
|
1538 |
|
subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, & |
1778 |
|
|
1779 |
|
function FF_RequiresPrepairCalc() result(doesit) |
1780 |
|
logical :: doesit |
1781 |
< |
doesit = FF_uses_EAM .or. FF_uses_SC & |
1732 |
< |
.or. FF_uses_MEAM |
1781 |
> |
doesit = FF_uses_EAM .or. FF_uses_SC |
1782 |
|
end function FF_RequiresPrepairCalc |
1783 |
|
|
1784 |
|
#ifdef PROFILE |