ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2266 by chuckv, Thu Jul 28 22:12:45 2005 UTC vs.
Revision 3131 by chrisfen, Wed May 2 00:18:08 2007 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
48	>	!! @version $Id: doForces.F90,v 1.88 2007-05-02 00:18:08 chrisfen Exp $, $Date: 2007-05-02 00:18:08 $, $Name: not supported by cvs2svn $, $Revision: 1.88 $
49
50
51		module doForces
#	Line 58 | Line 58 | module doForces
58		use lj
59		use sticky
60		use electrostatic_module
61	<	use reaction_field
62	<	use gb_pair
61	>	use gayberne
62		use shapes
63		use vector_class
64		use eam
65	+	use suttonchen
66		use status
67		#ifdef IS_MPI
68		use mpiSimulation
#	Line 72 | Line 72 | module doForces
72		PRIVATE
73
74		#define __FORTRAN90
75	<	#include "UseTheForce/fSwitchingFunction.h"
75	>	#include "UseTheForce/fCutoffPolicy.h"
76		#include "UseTheForce/DarkSide/fInteractionMap.h"
77	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78
79		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80		INTEGER, PARAMETER:: PAIR_LOOP    = 2
81
81	–	logical, save :: haveRlist = .false.
82		logical, save :: haveNeighborList = .false.
83		logical, save :: haveSIMvariables = .false.
84		logical, save :: haveSaneForceField = .false.
85	<	logical, save :: haveInteractionMap = .false.
85	>	logical, save :: haveInteractionHash = .false.
86	>	logical, save :: haveGtypeCutoffMap = .false.
87	>	logical, save :: haveDefaultCutoffs = .false.
88	>	logical, save :: haveSkinThickness = .false.
89	>	logical, save :: haveElectrostaticSummationMethod = .false.
90	>	logical, save :: haveCutoffPolicy = .false.
91	>	logical, save :: VisitCutoffsAfterComputing = .false.
92	>	logical, save :: do_box_dipole = .false.
93
94		logical, save :: FF_uses_DirectionalAtoms
88	–	logical, save :: FF_uses_LennardJones
89	–	logical, save :: FF_uses_Electrostatics
90	–	logical, save :: FF_uses_Charges
95		logical, save :: FF_uses_Dipoles
92	–	logical, save :: FF_uses_Quadrupoles
93	–	logical, save :: FF_uses_Sticky
94	–	logical, save :: FF_uses_StickyPower
96		logical, save :: FF_uses_GayBerne
97		logical, save :: FF_uses_EAM
98	<	logical, save :: FF_uses_Shapes
99	<	logical, save :: FF_uses_FLARB
100	<	logical, save :: FF_uses_RF
98	>	logical, save :: FF_uses_SC
99	>	logical, save :: FF_uses_MEAM
100	>
101
102		logical, save :: SIM_uses_DirectionalAtoms
102	–	logical, save :: SIM_uses_LennardJones
103	–	logical, save :: SIM_uses_Electrostatics
104	–	logical, save :: SIM_uses_Charges
105	–	logical, save :: SIM_uses_Dipoles
106	–	logical, save :: SIM_uses_Quadrupoles
107	–	logical, save :: SIM_uses_Sticky
108	–	logical, save :: SIM_uses_StickyPower
109	–	logical, save :: SIM_uses_GayBerne
103		logical, save :: SIM_uses_EAM
104	<	logical, save :: SIM_uses_Shapes
105	<	logical, save :: SIM_uses_FLARB
113	<	logical, save :: SIM_uses_RF
104	>	logical, save :: SIM_uses_SC
105	>	logical, save :: SIM_uses_MEAM
106		logical, save :: SIM_requires_postpair_calc
107		logical, save :: SIM_requires_prepair_calc
108		logical, save :: SIM_uses_PBC
109	<	logical, save :: SIM_uses_molecular_cutoffs
109	>	logical, save :: SIM_uses_AtomicVirial
110
111	<	!!!GO AWAY---------
112	<	!!!!!real(kind=dp), save :: rlist, rlistsq
111	>	integer, save :: electrostaticSummationMethod
112	>	integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
113
114	+	real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
115	+	real(kind=dp), save :: skinThickness
116	+	logical, save :: defaultDoShiftPot
117	+	logical, save :: defaultDoShiftFrc
118	+
119		public :: init_FF
120	+	public :: setCutoffs
121	+	public :: cWasLame
122	+	public :: setElectrostaticMethod
123	+	public :: setBoxDipole
124	+	public :: getBoxDipole
125	+	public :: setCutoffPolicy
126	+	public :: setSkinThickness
127		public :: do_force_loop
124	–	!  public :: setRlistDF
125	–	!public :: addInteraction
126	–	!public :: setInteractionHash
127	–	!public :: getInteractionHash
128	–	public :: createInteractionMap
129	–	public :: createRcuts
128
129		#ifdef PROFILE
130		public :: getforcetime
#	Line 134 | Line 132 | module doForces
132		real :: forceTimeInitial, forceTimeFinal
133		integer :: nLoops
134		#endif
137	–
138	–	type, public :: Interaction
139	–	integer :: InteractionHash
140	–	real(kind=dp) :: rList = 0.0_dp
141	–	real(kind=dp) :: rListSq = 0.0_dp
142	–	end type Interaction
135
136	<	type(Interaction), dimension(:,:),allocatable :: InteractionMap
137	<
136	>	!! Variables for cutoff mapping and interaction mapping
137	>	! Bit hash to determine pair-pair interactions.
138	>	integer, dimension(:,:), allocatable :: InteractionHash
139	>	real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
140	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
141	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
142
143	<
143	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
144	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
145	>
146	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
147	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
148	>	type ::gtypeCutoffs
149	>	real(kind=dp) :: rcut
150	>	real(kind=dp) :: rcutsq
151	>	real(kind=dp) :: rlistsq
152	>	end type gtypeCutoffs
153	>	type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
154	>
155	>	real(kind=dp), dimension(3) :: boxDipole
156	>
157		contains
158
159	<
151	<	subroutine createInteractionMap(status)
159	>	subroutine createInteractionHash()
160		integer :: nAtypes
153	–	integer, intent(out) :: status
161		integer :: i
162		integer :: j
163	<	integer :: ihash
164	<	real(kind=dp) :: myRcut
158	<	! Test Types
163	>	integer :: iHash
164	>	!! Test Types
165		logical :: i_is_LJ
166		logical :: i_is_Elect
167		logical :: i_is_Sticky
#	Line 163 | Line 169 | contains
169		logical :: i_is_GB
170		logical :: i_is_EAM
171		logical :: i_is_Shape
172	+	logical :: i_is_SC
173	+	logical :: i_is_MEAM
174		logical :: j_is_LJ
175		logical :: j_is_Elect
176		logical :: j_is_Sticky
#	Line 170 | Line 178 | contains
178		logical :: j_is_GB
179		logical :: j_is_EAM
180		logical :: j_is_Shape
181	<
182	<	status = 0
183	<
181	>	logical :: j_is_SC
182	>	logical :: j_is_MEAM
183	>	real(kind=dp) :: myRcut
184	>
185		if (.not. associated(atypes)) then
186	<	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
178	<	status = -1
186	>	call handleError("doForces", "atypes was not present before call of createInteractionHash!")
187		return
188		endif
189
190		nAtypes = getSize(atypes)
191
192		if (nAtypes == 0) then
193	<	status = -1
193	>	call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
194		return
195		end if
196
197	<	if (.not. allocated(InteractionMap)) then
198	<	allocate(InteractionMap(nAtypes,nAtypes))
197	>	if (.not. allocated(InteractionHash)) then
198	>	allocate(InteractionHash(nAtypes,nAtypes))
199	>	else
200	>	deallocate(InteractionHash)
201	>	allocate(InteractionHash(nAtypes,nAtypes))
202		endif
203	+
204	+	if (.not. allocated(atypeMaxCutoff)) then
205	+	allocate(atypeMaxCutoff(nAtypes))
206	+	else
207	+	deallocate(atypeMaxCutoff)
208	+	allocate(atypeMaxCutoff(nAtypes))
209	+	endif
210
211		do i = 1, nAtypes
212		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 198 | Line 216 | contains
216		call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
217		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
218		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
219	+	call getElementProperty(atypes, i, "is_SC", i_is_SC)
220	+	call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
221
222		do j = i, nAtypes
223
#	Line 211 | Line 231 | contains
231		call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
232		call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
233		call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
234	+	call getElementProperty(atypes, j, "is_SC", j_is_SC)
235	+	call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
236
237		if (i_is_LJ .and. j_is_LJ) then
238		iHash = ior(iHash, LJ_PAIR)
#	Line 232 | Line 254 | contains
254		iHash = ior(iHash, EAM_PAIR)
255		endif
256
257	+	if (i_is_SC .and. j_is_SC) then
258	+	iHash = ior(iHash, SC_PAIR)
259	+	endif
260	+
261		if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
262		if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
263		if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
#	Line 241 | Line 267 | contains
267		if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
268
269
270	<	InteractionMap(i,j)%InteractionHash = iHash
271	<	InteractionMap(j,i)%InteractionHash = iHash
270	>	InteractionHash(i,j) = iHash
271	>	InteractionHash(j,i) = iHash
272
273		end do
274
275		end do
250	–	end subroutine createInteractionMap
276
277	<	! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
278	<	subroutine createRcuts(defaultRList,stat)
254	<	real(kind=dp), intent(in), optional :: defaultRList
255	<	integer :: iMap
256	<	integer :: map_i,map_j
257	<	real(kind=dp) :: thisRCut = 0.0_dp
258	<	real(kind=dp) :: actualCutoff = 0.0_dp
259	<	integer, intent(out) :: stat
260	<	integer :: nAtypes
261	<	integer :: myStatus
277	>	haveInteractionHash = .true.
278	>	end subroutine createInteractionHash
279
280	<	stat = 0
264	<	if (.not. haveInteractionMap) then
280	>	subroutine createGtypeCutoffMap()
281
282	<	call createInteractionMap(myStatus)
282	>	logical :: i_is_LJ
283	>	logical :: i_is_Elect
284	>	logical :: i_is_Sticky
285	>	logical :: i_is_StickyP
286	>	logical :: i_is_GB
287	>	logical :: i_is_EAM
288	>	logical :: i_is_Shape
289	>	logical :: i_is_SC
290	>	logical :: GtypeFound
291
292	<	if (myStatus .ne. 0) then
293	<	write(default_error, *) 'createInteractionMap failed in doForces!'
294	<	stat = -1
295	<	return
292	>	integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
293	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
294	>	integer :: nGroupsInRow
295	>	integer :: nGroupsInCol
296	>	integer :: nGroupTypesRow,nGroupTypesCol
297	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
298	>	real(kind=dp) :: biggestAtypeCutoff
299	>
300	>	if (.not. haveInteractionHash) then
301	>	call createInteractionHash()
302	>	endif
303	>	#ifdef IS_MPI
304	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
305	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
306	>	#endif
307	>	nAtypes = getSize(atypes)
308	>	! Set all of the initial cutoffs to zero.
309	>	atypeMaxCutoff = 0.0_dp
310	>	do i = 1, nAtypes
311	>	if (SimHasAtype(i)) then
312	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
313	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
314	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
315	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
316	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
317	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
318	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
319	>	call getElementProperty(atypes, i, "is_SC", i_is_SC)
320	>
321	>	if (haveDefaultCutoffs) then
322	>	atypeMaxCutoff(i) = defaultRcut
323	>	else
324	>	if (i_is_LJ) then
325	>	thisRcut = getSigma(i) * 2.5_dp
326	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327	>	endif
328	>	if (i_is_Elect) then
329	>	thisRcut = defaultRcut
330	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331	>	endif
332	>	if (i_is_Sticky) then
333	>	thisRcut = getStickyCut(i)
334	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
335	>	endif
336	>	if (i_is_StickyP) then
337	>	thisRcut = getStickyPowerCut(i)
338	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
339	>	endif
340	>	if (i_is_GB) then
341	>	thisRcut = getGayBerneCut(i)
342	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
343	>	endif
344	>	if (i_is_EAM) then
345	>	thisRcut = getEAMCut(i)
346	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
347	>	endif
348	>	if (i_is_Shape) then
349	>	thisRcut = getShapeCut(i)
350	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
351	>	endif
352	>	if (i_is_SC) then
353	>	thisRcut = getSCCut(i)
354	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
355	>	endif
356	>	endif
357	>
358	>	if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
359	>	biggestAtypeCutoff = atypeMaxCutoff(i)
360	>	endif
361	>
362		endif
363	+	enddo
364	+
365	+	istart = 1
366	+	jstart = 1
367	+	#ifdef IS_MPI
368	+	iend = nGroupsInRow
369	+	jend = nGroupsInCol
370	+	#else
371	+	iend = nGroups
372	+	jend = nGroups
373	+	#endif
374	+
375	+	!! allocate the groupToGtype and gtypeMaxCutoff here.
376	+	if(.not.allocated(groupToGtypeRow)) then
377	+	!  allocate(groupToGtype(iend))
378	+	allocate(groupToGtypeRow(iend))
379	+	else
380	+	deallocate(groupToGtypeRow)
381	+	allocate(groupToGtypeRow(iend))
382		endif
383	+	if(.not.allocated(groupMaxCutoffRow)) then
384	+	allocate(groupMaxCutoffRow(iend))
385	+	else
386	+	deallocate(groupMaxCutoffRow)
387	+	allocate(groupMaxCutoffRow(iend))
388	+	end if
389
390	+	if(.not.allocated(gtypeMaxCutoffRow)) then
391	+	allocate(gtypeMaxCutoffRow(iend))
392	+	else
393	+	deallocate(gtypeMaxCutoffRow)
394	+	allocate(gtypeMaxCutoffRow(iend))
395	+	endif
396
397	<	nAtypes = getSize(atypes)
398	<	! If we pass a default rcut, set all atypes to that cutoff distance
399	<	if(present(defaultRList)) then
400	<	InteractionMap(:,:)%rList = defaultRList
401	<	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
402	<	haveRlist = .true.
403	<	return
397	>
398	>	#ifdef IS_MPI
399	>	! We only allocate new storage if we are in MPI because Ncol /= Nrow
400	>	if(.not.associated(groupToGtypeCol)) then
401	>	allocate(groupToGtypeCol(jend))
402	>	else
403	>	deallocate(groupToGtypeCol)
404	>	allocate(groupToGtypeCol(jend))
405		end if
406
407	<	do map_i = 1,nAtypes
408	<	do map_j = map_i,nAtypes
409	<	iMap = InteractionMap(map_i, map_j)%InteractionHash
410	<
411	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
412	<	!            thisRCut = getLJCutOff(map_i,map_j)
413	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
414	<	endif
415	<
416	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
417	<	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
418	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
297	<	endif
298	<
299	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300	<	!             thisRCut = getStickyCutOff(map_i,map_j)
301	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
302	<	endif
303	<
304	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305	<	!              thisRCut = getStickyPowerCutOff(map_i,map_j)
306	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
307	<	endif
308	<
309	<	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310	<	!              thisRCut = getGayberneCutOff(map_i,map_j)
311	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
312	<	endif
313	<
314	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315	<	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
317	<	endif
318	<
319	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
320	<	!              thisRCut = getEAMCutOff(map_i,map_j)
321	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
322	<	endif
323	<
324	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
325	<	!              thisRCut = getShapeCutOff(map_i,map_j)
326	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
327	<	endif
328	<
329	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
330	<	!              thisRCut = getShapeLJCutOff(map_i,map_j)
331	<	if (thisRcut > actualCutoff) actualCutoff = thisRcut
332	<	endif
333	<	InteractionMap(map_i, map_j)%rList = actualCutoff
334	<	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335	<	end do
336	<	end do
337	<	haveRlist = .true.
338	<	end subroutine createRcuts
407	>	if(.not.associated(groupMaxCutoffCol)) then
408	>	allocate(groupMaxCutoffCol(jend))
409	>	else
410	>	deallocate(groupMaxCutoffCol)
411	>	allocate(groupMaxCutoffCol(jend))
412	>	end if
413	>	if(.not.associated(gtypeMaxCutoffCol)) then
414	>	allocate(gtypeMaxCutoffCol(jend))
415	>	else
416	>	deallocate(gtypeMaxCutoffCol)
417	>	allocate(gtypeMaxCutoffCol(jend))
418	>	end if
419
420	+	groupMaxCutoffCol = 0.0_dp
421	+	gtypeMaxCutoffCol = 0.0_dp
422
423	<	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
424	<	!!$  subroutine setRlistDF( this_rlist )
425	<	!!$
344	<	!!$   real(kind=dp) :: this_rlist
345	<	!!$
346	<	!!$    rlist = this_rlist
347	<	!!$    rlistsq = rlist * rlist
348	<	!!$
349	<	!!$    haveRlist = .true.
350	<	!!$
351	<	!!$  end subroutine setRlistDF
423	>	#endif
424	>	groupMaxCutoffRow = 0.0_dp
425	>	gtypeMaxCutoffRow = 0.0_dp
426
427
428	<	subroutine setSimVariables()
429	<	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
430	<	SIM_uses_LennardJones = SimUsesLennardJones()
431	<	SIM_uses_Electrostatics = SimUsesElectrostatics()
432	<	SIM_uses_Charges = SimUsesCharges()
433	<	SIM_uses_Dipoles = SimUsesDipoles()
434	<	SIM_uses_Sticky = SimUsesSticky()
435	<	SIM_uses_StickyPower = SimUsesStickyPower()
436	<	SIM_uses_GayBerne = SimUsesGayBerne()
437	<	SIM_uses_EAM = SimUsesEAM()
438	<	SIM_uses_Shapes = SimUsesShapes()
439	<	SIM_uses_FLARB = SimUsesFLARB()
440	<	SIM_uses_RF = SimUsesRF()
441	<	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
442	<	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
443	<	SIM_uses_PBC = SimUsesPBC()
444	<
445	<	haveSIMvariables = .true.
428	>	!! first we do a single loop over the cutoff groups to find the
429	>	!! largest cutoff for any atypes present in this group.  We also
430	>	!! create gtypes at this point.
431	>
432	>	tol = 1.0e-6_dp
433	>	nGroupTypesRow = 0
434	>	nGroupTypesCol = 0
435	>	do i = istart, iend
436	>	n_in_i = groupStartRow(i+1) - groupStartRow(i)
437	>	groupMaxCutoffRow(i) = 0.0_dp
438	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
439	>	atom1 = groupListRow(ia)
440	>	#ifdef IS_MPI
441	>	me_i = atid_row(atom1)
442	>	#else
443	>	me_i = atid(atom1)
444	>	#endif
445	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
446	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
447	>	endif
448	>	enddo
449	>	if (nGroupTypesRow.eq.0) then
450	>	nGroupTypesRow = nGroupTypesRow + 1
451	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
452	>	groupToGtypeRow(i) = nGroupTypesRow
453	>	else
454	>	GtypeFound = .false.
455	>	do g = 1, nGroupTypesRow
456	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
457	>	groupToGtypeRow(i) = g
458	>	GtypeFound = .true.
459	>	endif
460	>	enddo
461	>	if (.not.GtypeFound) then
462	>	nGroupTypesRow = nGroupTypesRow + 1
463	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
464	>	groupToGtypeRow(i) = nGroupTypesRow
465	>	endif
466	>	endif
467	>	enddo
468
469	<	return
470	<	end subroutine setSimVariables
469	>	#ifdef IS_MPI
470	>	do j = jstart, jend
471	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
472	>	groupMaxCutoffCol(j) = 0.0_dp
473	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
474	>	atom1 = groupListCol(ja)
475
476	<	subroutine doReadyCheck(error)
377	<	integer, intent(out) :: error
476	>	me_j = atid_col(atom1)
477
478	<	integer :: myStatus
478	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
479	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
480	>	endif
481	>	enddo
482
483	<	error = 0
483	>	if (nGroupTypesCol.eq.0) then
484	>	nGroupTypesCol = nGroupTypesCol + 1
485	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
486	>	groupToGtypeCol(j) = nGroupTypesCol
487	>	else
488	>	GtypeFound = .false.
489	>	do g = 1, nGroupTypesCol
490	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
491	>	groupToGtypeCol(j) = g
492	>	GtypeFound = .true.
493	>	endif
494	>	enddo
495	>	if (.not.GtypeFound) then
496	>	nGroupTypesCol = nGroupTypesCol + 1
497	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
498	>	groupToGtypeCol(j) = nGroupTypesCol
499	>	endif
500	>	endif
501	>	enddo
502
503	<	if (.not. haveInteractionMap) then
503	>	#else
504	>	! Set pointers to information we just found
505	>	nGroupTypesCol = nGroupTypesRow
506	>	groupToGtypeCol => groupToGtypeRow
507	>	gtypeMaxCutoffCol => gtypeMaxCutoffRow
508	>	groupMaxCutoffCol => groupMaxCutoffRow
509	>	#endif
510
511	<	myStatus = 0
511	>	!! allocate the gtypeCutoffMap here.
512	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
513	>	!! then we do a double loop over all the group TYPES to find the cutoff
514	>	!! map between groups of two types
515	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
516
517	<	call createInteractionMap(myStatus)
517	>	do i = 1, nGroupTypesRow
518	>	do j = 1, nGroupTypesCol
519	>
520	>	select case(cutoffPolicy)
521	>	case(TRADITIONAL_CUTOFF_POLICY)
522	>	thisRcut = tradRcut
523	>	case(MIX_CUTOFF_POLICY)
524	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
525	>	case(MAX_CUTOFF_POLICY)
526	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
527	>	case default
528	>	call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
529	>	return
530	>	end select
531	>	gtypeCutoffMap(i,j)%rcut = thisRcut
532	>
533	>	if (thisRcut.gt.largestRcut) largestRcut = thisRcut
534
535	<	if (myStatus .ne. 0) then
536	<	write(default_error, *) 'createInteractionMap failed in doForces!'
537	<	error = -1
538	<	return
539	<	endif
535	>	gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
536	>
537	>	if (.not.haveSkinThickness) then
538	>	skinThickness = 1.0_dp
539	>	endif
540	>
541	>	gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
542	>
543	>	! sanity check
544	>
545	>	if (haveDefaultCutoffs) then
546	>	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
547	>	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
548	>	endif
549	>	endif
550	>	enddo
551	>	enddo
552	>
553	>	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
554	>	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
555	>	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
556	>	#ifdef IS_MPI
557	>	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
558	>	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
559	>	#endif
560	>	groupMaxCutoffCol => null()
561	>	gtypeMaxCutoffCol => null()
562	>
563	>	haveGtypeCutoffMap = .true.
564	>	end subroutine createGtypeCutoffMap
565	>
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	>	defaultDoShiftPot = defSP
577	>	defaultDoShiftFrc = defSF
578	>
579	>	if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
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, defaultDoShiftPot, &
602	>	defaultDoShiftFrc )
603	>	call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
604	>	call setCutoffEAM( defaultRcut )
605	>	call setCutoffSC( defaultRcut )
606	>	call set_switch(defaultRsw, defaultRcut)
607	>	call setHmatDangerousRcutValue(defaultRcut)
608	>
609	>	haveDefaultCutoffs = .true.
610	>	haveGtypeCutoffMap = .false.
611	>
612	>	end subroutine setCutoffs
613	>
614	>	subroutine cWasLame()
615	>
616	>	VisitCutoffsAfterComputing = .true.
617	>	return
618	>
619	>	end subroutine cWasLame
620	>
621	>	subroutine setCutoffPolicy(cutPolicy)
622	>
623	>	integer, intent(in) :: cutPolicy
624	>
625	>	cutoffPolicy = cutPolicy
626	>	haveCutoffPolicy = .true.
627	>	haveGtypeCutoffMap = .false.
628	>
629	>	end subroutine setCutoffPolicy
630	>
631	>	subroutine setBoxDipole()
632	>
633	>	do_box_dipole = .true.
634	>
635	>	end subroutine setBoxDipole
636	>
637	>	subroutine getBoxDipole( box_dipole )
638	>
639	>	real(kind=dp), intent(inout), dimension(3) :: box_dipole
640	>
641	>	box_dipole = boxDipole
642	>
643	>	end subroutine getBoxDipole
644	>
645	>	subroutine setElectrostaticMethod( thisESM )
646	>
647	>	integer, intent(in) :: thisESM
648	>
649	>	electrostaticSummationMethod = thisESM
650	>	haveElectrostaticSummationMethod = .true.
651	>
652	>	end subroutine setElectrostaticMethod
653	>
654	>	subroutine setSkinThickness( thisSkin )
655	>
656	>	real(kind=dp), intent(in) :: thisSkin
657	>
658	>	skinThickness = thisSkin
659	>	haveSkinThickness = .true.
660	>	haveGtypeCutoffMap = .false.
661	>
662	>	end subroutine setSkinThickness
663	>
664	>	subroutine setSimVariables()
665	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
666	>	SIM_uses_EAM = SimUsesEAM()
667	>	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
668	>	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
669	>	SIM_uses_PBC = SimUsesPBC()
670	>	SIM_uses_SC = SimUsesSC()
671	>	SIM_uses_AtomicVirial = SimUsesAtomicVirial()
672	>
673	>	haveSIMvariables = .true.
674	>
675	>	return
676	>	end subroutine setSimVariables
677	>
678	>	subroutine doReadyCheck(error)
679	>	integer, intent(out) :: error
680	>	integer :: myStatus
681	>
682	>	error = 0
683	>
684	>	if (.not. haveInteractionHash) then
685	>	call createInteractionHash()
686		endif
687
688	<	if (.not. haveSIMvariables) then
689	<	call setSimVariables()
688	>	if (.not. haveGtypeCutoffMap) then
689	>	call createGtypeCutoffMap()
690		endif
691
692	<	if (.not. haveRlist) then
693	<	write(default_error, *) 'rList has not been set in doForces!'
694	<	error = -1
695	<	return
692	>	if (VisitCutoffsAfterComputing) then
693	>	call set_switch(largestRcut, largestRcut)
694	>	call setHmatDangerousRcutValue(largestRcut)
695	>	call setCutoffEAM(largestRcut)
696	>	call setCutoffSC(largestRcut)
697	>	VisitCutoffsAfterComputing = .false.
698		endif
699
700	+	if (.not. haveSIMvariables) then
701	+	call setSimVariables()
702	+	endif
703	+
704		if (.not. haveNeighborList) then
705		write(default_error, *) 'neighbor list has not been initialized in doForces!'
706		error = -1
707		return
708		end if
709	<
709	>
710		if (.not. haveSaneForceField) then
711		write(default_error, *) 'Force Field is not sane in doForces!'
712		error = -1
713		return
714		end if
715	<
715	>
716		#ifdef IS_MPI
717		if (.not. isMPISimSet()) then
718		write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
#	Line 426 | Line 724 | contains
724		end subroutine doReadyCheck
725
726
727	<	subroutine init_FF(use_RF_c, thisStat)
727	>	subroutine init_FF(thisStat)
728
431	–	logical, intent(in) :: use_RF_c
432	–
729		integer, intent(out) :: thisStat
730		integer :: my_status, nMatches
731		integer, pointer :: MatchList(:) => null()
436	–	real(kind=dp) :: rcut, rrf, rt, dielect
732
733		!! assume things are copacetic, unless they aren't
734		thisStat = 0
735
441	–	!! Fortran's version of a cast:
442	–	FF_uses_RF = use_RF_c
443	–
736		!! init_FF is called *after* all of the atom types have been
737		!! defined in atype_module using the new_atype subroutine.
738		!!
#	Line 448 | Line 740 | contains
740		!! interactions are used by the force field.
741
742		FF_uses_DirectionalAtoms = .false.
451	–	FF_uses_LennardJones = .false.
452	–	FF_uses_Electrostatics = .false.
453	–	FF_uses_Charges = .false.
743		FF_uses_Dipoles = .false.
455	–	FF_uses_Sticky = .false.
456	–	FF_uses_StickyPower = .false.
744		FF_uses_GayBerne = .false.
745		FF_uses_EAM = .false.
746	<	FF_uses_Shapes = .false.
460	<	FF_uses_FLARB = .false.
746	>	FF_uses_SC = .false.
747
748		call getMatchingElementList(atypes, "is_Directional", .true., &
749		nMatches, MatchList)
750		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
751
466	–	call getMatchingElementList(atypes, "is_LennardJones", .true., &
467	–	nMatches, MatchList)
468	–	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
469	–
470	–	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
471	–	nMatches, MatchList)
472	–	if (nMatches .gt. 0) then
473	–	FF_uses_Electrostatics = .true.
474	–	endif
475	–
476	–	call getMatchingElementList(atypes, "is_Charge", .true., &
477	–	nMatches, MatchList)
478	–	if (nMatches .gt. 0) then
479	–	FF_uses_Charges = .true.
480	–	FF_uses_Electrostatics = .true.
481	–	endif
482	–
752		call getMatchingElementList(atypes, "is_Dipole", .true., &
753		nMatches, MatchList)
754	<	if (nMatches .gt. 0) then
486	<	FF_uses_Dipoles = .true.
487	<	FF_uses_Electrostatics = .true.
488	<	FF_uses_DirectionalAtoms = .true.
489	<	endif
490	<
491	<	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
492	<	nMatches, MatchList)
493	<	if (nMatches .gt. 0) then
494	<	FF_uses_Quadrupoles = .true.
495	<	FF_uses_Electrostatics = .true.
496	<	FF_uses_DirectionalAtoms = .true.
497	<	endif
498	<
499	<	call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
500	<	MatchList)
501	<	if (nMatches .gt. 0) then
502	<	FF_uses_Sticky = .true.
503	<	FF_uses_DirectionalAtoms = .true.
504	<	endif
505	<
506	<	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
507	<	MatchList)
508	<	if (nMatches .gt. 0) then
509	<	FF_uses_StickyPower = .true.
510	<	FF_uses_DirectionalAtoms = .true.
511	<	endif
754	>	if (nMatches .gt. 0) FF_uses_Dipoles = .true.
755
756		call getMatchingElementList(atypes, "is_GayBerne", .true., &
757		nMatches, MatchList)
758	<	if (nMatches .gt. 0) then
516	<	FF_uses_GayBerne = .true.
517	<	FF_uses_DirectionalAtoms = .true.
518	<	endif
758	>	if (nMatches .gt. 0) FF_uses_GayBerne = .true.
759
760		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
761		if (nMatches .gt. 0) FF_uses_EAM = .true.
762
763	<	call getMatchingElementList(atypes, "is_Shape", .true., &
764	<	nMatches, MatchList)
525	<	if (nMatches .gt. 0) then
526	<	FF_uses_Shapes = .true.
527	<	FF_uses_DirectionalAtoms = .true.
528	<	endif
763	>	call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
764	>	if (nMatches .gt. 0) FF_uses_SC = .true.
765
530	–	call getMatchingElementList(atypes, "is_FLARB", .true., &
531	–	nMatches, MatchList)
532	–	if (nMatches .gt. 0) FF_uses_FLARB = .true.
766
534	–	!! Assume sanity (for the sake of argument)
767		haveSaneForceField = .true.
768
537	–	!! check to make sure the FF_uses_RF setting makes sense
538	–
539	–	if (FF_uses_dipoles) then
540	–	if (FF_uses_RF) then
541	–	dielect = getDielect()
542	–	call initialize_rf(dielect)
543	–	endif
544	–	else
545	–	if (FF_uses_RF) then
546	–	write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
547	–	thisStat = -1
548	–	haveSaneForceField = .false.
549	–	return
550	–	endif
551	–	endif
552	–
553	–	!sticky module does not contain check_sticky_FF anymore
554	–	!if (FF_uses_sticky) then
555	–	!   call check_sticky_FF(my_status)
556	–	!   if (my_status /= 0) then
557	–	!      thisStat = -1
558	–	!      haveSaneForceField = .false.
559	–	!      return
560	–	!   end if
561	–	!endif
562	–
769		if (FF_uses_EAM) then
770		call init_EAM_FF(my_status)
771		if (my_status /= 0) then
#	Line 570 | Line 776 | contains
776		end if
777		endif
778
573	–	if (FF_uses_GayBerne) then
574	–	call check_gb_pair_FF(my_status)
575	–	if (my_status .ne. 0) then
576	–	thisStat = -1
577	–	haveSaneForceField = .false.
578	–	return
579	–	endif
580	–	endif
581	–
582	–	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583	–	endif
584	–
779		if (.not. haveNeighborList) then
780		!! Create neighbor lists
781		call expandNeighborList(nLocal, my_status)
#	Line 615 | Line 809 | contains
809
810		!! Stress Tensor
811		real( kind = dp), dimension(9) :: tau
812	<	real ( kind = dp ) :: pot
812	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
813		logical ( kind = 2) :: do_pot_c, do_stress_c
814		logical :: do_pot
815		logical :: do_stress
816		logical :: in_switching_region
817		#ifdef IS_MPI
818	<	real( kind = DP ) :: pot_local
818	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
819		integer :: nAtomsInRow
820		integer :: nAtomsInCol
821		integer :: nprocs
#	Line 634 | Line 828 | contains
828		integer :: istart, iend
829		integer :: ia, jb, atom1, atom2
830		integer :: nlist
831	<	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
831	>	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
832		real( kind = DP ) :: sw, dswdr, swderiv, mf
833	<	real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
834	<	real(kind=dp) :: rfpot, mu_i, virial
833	>	real( kind = DP ) :: rVal
834	>	real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
835	>	real(kind=dp) :: rfpot, mu_i
836	>	real(kind=dp):: rCut
837		integer :: me_i, me_j, n_in_i, n_in_j
838		logical :: is_dp_i
839		integer :: neighborListSize
#	Line 645 | Line 841 | contains
841		integer :: localError
842		integer :: propPack_i, propPack_j
843		integer :: loopStart, loopEnd, loop
844	<	integer :: iMap
845	<	real(kind=dp) :: listSkin = 1.0
844	>	integer :: iHash
845	>	integer :: i1
846
847	+	!! the variables for the box dipole moment
848	+	#ifdef IS_MPI
849	+	integer :: pChgCount_local
850	+	integer :: nChgCount_local
851	+	real(kind=dp) :: pChg_local
852	+	real(kind=dp) :: nChg_local
853	+	real(kind=dp), dimension(3) :: pChgPos_local
854	+	real(kind=dp), dimension(3) :: nChgPos_local
855	+	real(kind=dp), dimension(3) :: dipVec_local
856	+	#endif
857	+	integer :: pChgCount
858	+	integer :: nChgCount
859	+	real(kind=dp) :: pChg
860	+	real(kind=dp) :: nChg
861	+	real(kind=dp) :: chg_value
862	+	real(kind=dp), dimension(3) :: pChgPos
863	+	real(kind=dp), dimension(3) :: nChgPos
864	+	real(kind=dp), dimension(3) :: dipVec
865	+	real(kind=dp), dimension(3) :: chgVec
866	+
867	+	!! initialize box dipole variables
868	+	if (do_box_dipole) then
869	+	#ifdef IS_MPI
870	+	pChg_local = 0.0_dp
871	+	nChg_local = 0.0_dp
872	+	pChgCount_local = 0
873	+	nChgCount_local = 0
874	+	do i=1, 3
875	+	pChgPos_local = 0.0_dp
876	+	nChgPos_local = 0.0_dp
877	+	dipVec_local = 0.0_dp
878	+	enddo
879	+	#endif
880	+	pChg = 0.0_dp
881	+	nChg = 0.0_dp
882	+	pChgCount = 0
883	+	nChgCount = 0
884	+	chg_value = 0.0_dp
885	+
886	+	do i=1, 3
887	+	pChgPos(i) = 0.0_dp
888	+	nChgPos(i) = 0.0_dp
889	+	dipVec(i) = 0.0_dp
890	+	chgVec(i) = 0.0_dp
891	+	boxDipole(i) = 0.0_dp
892	+	enddo
893	+	endif
894	+
895		!! initialize local variables
896
897		#ifdef IS_MPI
#	Line 710 | Line 954 | contains
954		! (but only on the first time through):
955		if (loop .eq. loopStart) then
956		#ifdef IS_MPI
957	<	call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
957	>	call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
958		update_nlist)
959		#else
960	<	call checkNeighborList(nGroups, q_group, listSkin, &
960	>	call checkNeighborList(nGroups, q_group, skinThickness, &
961		update_nlist)
962		#endif
963		endif
#	Line 743 | Line 987 | contains
987
988		if (update_nlist) then
989		#ifdef IS_MPI
746	–	me_i = atid_row(i)
990		jstart = 1
991		jend = nGroupsInCol
992		#else
750	–	me_i = atid(i)
993		jstart = i+1
994		jend = nGroups
995		#endif
#	Line 773 | Line 1015 | contains
1015		me_j = atid(j)
1016		call get_interatomic_vector(q_group(:,i), &
1017		q_group(:,j), d_grp, rgrpsq)
1018	<	#endif
1018	>	#endif
1019
1020	<	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
1020	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1021		if (update_nlist) then
1022		nlist = nlist + 1
1023
#	Line 795 | Line 1037 | contains
1037
1038		list(nlist) = j
1039		endif
1040	<
1041	<	if (loop .eq. PAIR_LOOP) then
800	<	vij = 0.0d0
801	<	fij(1:3) = 0.0d0
802	<	endif
803	<
804	<	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
805	<	in_switching_region)
806	<
807	<	n_in_j = groupStartCol(j+1) - groupStartCol(j)
808	<
809	<	do ia = groupStartRow(i), groupStartRow(i+1)-1
810	<
811	<	atom1 = groupListRow(ia)
812	<
813	<	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
814	<
815	<	atom2 = groupListCol(jb)
816	<
817	<	if (skipThisPair(atom1, atom2)) cycle inner
818	<
819	<	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
820	<	d_atm(1:3) = d_grp(1:3)
821	<	ratmsq = rgrpsq
822	<	else
823	<	#ifdef IS_MPI
824	<	call get_interatomic_vector(q_Row(:,atom1), &
825	<	q_Col(:,atom2), d_atm, ratmsq)
826	<	#else
827	<	call get_interatomic_vector(q(:,atom1), &
828	<	q(:,atom2), d_atm, ratmsq)
829	<	#endif
830	<	endif
1040	>
1041	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1042
1043	<	if (loop .eq. PREPAIR_LOOP) then
1043	>	rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1044	>	if (loop .eq. PAIR_LOOP) then
1045	>	vij = 0.0_dp
1046	>	fij(1) = 0.0_dp
1047	>	fij(2) = 0.0_dp
1048	>	fij(3) = 0.0_dp
1049	>	endif
1050	>
1051	>	call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1052	>
1053	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
1054	>
1055	>	do ia = groupStartRow(i), groupStartRow(i+1)-1
1056	>
1057	>	atom1 = groupListRow(ia)
1058	>
1059	>	inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1060	>
1061	>	atom2 = groupListCol(jb)
1062	>
1063	>	if (skipThisPair(atom1, atom2))  cycle inner
1064	>
1065	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1066	>	d_atm(1) = d_grp(1)
1067	>	d_atm(2) = d_grp(2)
1068	>	d_atm(3) = d_grp(3)
1069	>	ratmsq = rgrpsq
1070	>	else
1071	>	#ifdef IS_MPI
1072	>	call get_interatomic_vector(q_Row(:,atom1), &
1073	>	q_Col(:,atom2), d_atm, ratmsq)
1074	>	#else
1075	>	call get_interatomic_vector(q(:,atom1), &
1076	>	q(:,atom2), d_atm, ratmsq)
1077	>	#endif
1078	>	endif
1079	>
1080	>	if (loop .eq. PREPAIR_LOOP) then
1081		#ifdef IS_MPI
1082	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1083	<	rgrpsq, d_grp, do_pot, do_stress, &
1084	<	eFrame, A, f, t, pot_local)
1082	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1083	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1084	>	eFrame, A, f, t, pot_local)
1085		#else
1086	<	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1087	<	rgrpsq, d_grp, do_pot, do_stress, &
1088	<	eFrame, A, f, t, pot)
1086	>	call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1087	>	rgrpsq, d_grp, rCut, do_pot, do_stress, &
1088	>	eFrame, A, f, t, pot)
1089		#endif
1090	<	else
1090	>	else
1091		#ifdef IS_MPI
1092	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1093	<	do_pot, &
1094	<	eFrame, A, f, t, pot_local, vpair, fpair)
1092	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1093	>	do_pot, eFrame, A, f, t, pot_local, vpair, &
1094	>	fpair, d_grp, rgrp, rCut)
1095		#else
1096	<	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1097	<	do_pot,  &
1098	<	eFrame, A, f, t, pot, vpair, fpair)
1096	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1097	>	do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1098	>	d_grp, rgrp, rCut)
1099		#endif
1100	+	vij = vij + vpair
1101	+	fij(1) = fij(1) + fpair(1)
1102	+	fij(2) = fij(2) + fpair(2)
1103	+	fij(3) = fij(3) + fpair(3)
1104	+	if (do_stress) then
1105	+	call add_stress_tensor(d_atm, fpair, tau)
1106	+	endif
1107	+	endif
1108	+	enddo inner
1109	+	enddo
1110
1111	<	vij = vij + vpair
1112	<	fij(1:3) = fij(1:3) + fpair(1:3)
1113	<	endif
1114	<	enddo inner
1115	<	enddo
1116	<
1117	<	if (loop .eq. PAIR_LOOP) then
1118	<	if (in_switching_region) then
1119	<	swderiv = vij*dswdr/rgrp
1120	<	fij(1) = fij(1) + swderiv*d_grp(1)
1121	<	fij(2) = fij(2) + swderiv*d_grp(2)
1122	<	fij(3) = fij(3) + swderiv*d_grp(3)
1123	<
1124	<	do ia=groupStartRow(i), groupStartRow(i+1)-1
1125	<	atom1=groupListRow(ia)
1126	<	mf = mfactRow(atom1)
1111	>	if (loop .eq. PAIR_LOOP) then
1112	>	if (in_switching_region) then
1113	>	swderiv = vij*dswdr/rgrp
1114	>	fg = swderiv*d_grp
1115	>
1116	>	fij(1) = fij(1) + fg(1)
1117	>	fij(2) = fij(2) + fg(2)
1118	>	fij(3) = fij(3) + fg(3)
1119	>
1120	>	if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1121	>	call add_stress_tensor(d_atm, fg, tau)
1122	>	endif
1123	>
1124	>	do ia=groupStartRow(i), groupStartRow(i+1)-1
1125	>	atom1=groupListRow(ia)
1126	>	mf = mfactRow(atom1)
1127	>	! fg is the force on atom ia due to cutoff group's
1128	>	! presence in switching region
1129	>	fg = swderiv*d_grp*mf
1130		#ifdef IS_MPI
1131	<	f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1132	<	f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1133	<	f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1131	>	f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1132	>	f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1133	>	f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1134		#else
1135	<	f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1136	<	f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1137	<	f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1135	>	f(1,atom1) = f(1,atom1) + fg(1)
1136	>	f(2,atom1) = f(2,atom1) + fg(2)
1137	>	f(3,atom1) = f(3,atom1) + fg(3)
1138		#endif
1139	<	enddo
1140	<
1141	<	do jb=groupStartCol(j), groupStartCol(j+1)-1
1142	<	atom2=groupListCol(jb)
882	<	mf = mfactCol(atom2)
1139	>	if (n_in_i .gt. 1) then
1140	>	if (do_stress.and.SIM_uses_AtomicVirial) then
1141	>	! find the distance between the atom and the center of
1142	>	! the cutoff group:
1143		#ifdef IS_MPI
1144	<	f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1145	<	f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
886	<	f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1144	>	call get_interatomic_vector(q_Row(:,atom1), &
1145	>	q_group_Row(:,i), dag, rag)
1146		#else
1147	<	f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1148	<	f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
890	<	f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1147	>	call get_interatomic_vector(q(:,atom1), &
1148	>	q_group(:,i), dag, rag)
1149		#endif
1150	<	enddo
1150	>	call add_stress_tensor(dag,fg,tau)
1151	>	endif
1152	>	endif
1153	>	enddo
1154	>
1155	>	do jb=groupStartCol(j), groupStartCol(j+1)-1
1156	>	atom2=groupListCol(jb)
1157	>	mf = mfactCol(atom2)
1158	>	! fg is the force on atom jb due to cutoff group's
1159	>	! presence in switching region
1160	>	fg = -swderiv*d_grp*mf
1161	>	#ifdef IS_MPI
1162	>	f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1163	>	f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1164	>	f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1165	>	#else
1166	>	f(1,atom2) = f(1,atom2) + fg(1)
1167	>	f(2,atom2) = f(2,atom2) + fg(2)
1168	>	f(3,atom2) = f(3,atom2) + fg(3)
1169	>	#endif
1170	>	if (n_in_j .gt. 1) then
1171	>	if (do_stress.and.SIM_uses_AtomicVirial) then
1172	>	! find the distance between the atom and the center of
1173	>	! the cutoff group:
1174	>	#ifdef IS_MPI
1175	>	call get_interatomic_vector(q_Col(:,atom2), &
1176	>	q_group_Col(:,j), dag, rag)
1177	>	#else
1178	>	call get_interatomic_vector(q(:,atom2), &
1179	>	q_group(:,j), dag, rag)
1180	>	#endif
1181	>	call add_stress_tensor(dag,fg,tau)
1182	>	endif
1183	>	endif
1184	>	enddo
1185	>	endif
1186		endif
894	–
895	–	if (do_stress) call add_stress_tensor(d_grp, fij)
1187		endif
1188	<	end if
1188	>	endif
1189		enddo
1190	+
1191		enddo outer
1192
1193		if (update_nlist) then
#	Line 955 | Line 1247 | contains
1247
1248		if (do_pot) then
1249		! scatter/gather pot_row into the members of my column
1250	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1251	<
1250	>	do i = 1,LR_POT_TYPES
1251	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1252	>	end do
1253		! scatter/gather pot_local into all other procs
1254		! add resultant to get total pot
1255		do i = 1, nlocal
1256	<	pot_local = pot_local + pot_Temp(i)
1256	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1257	>	+ pot_Temp(1:LR_POT_TYPES,i)
1258		enddo
1259
1260		pot_Temp = 0.0_DP
1261	<
1262	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1261	>	do i = 1,LR_POT_TYPES
1262	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1263	>	end do
1264		do i = 1, nlocal
1265	<	pot_local = pot_local + pot_Temp(i)
1265	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1266	>	+ pot_Temp(1:LR_POT_TYPES,i)
1267		enddo
1268
1269		endif
1270		#endif
1271
1272	<	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1272	>	if (SIM_requires_postpair_calc) then
1273	>	do i = 1, nlocal
1274	>
1275	>	! we loop only over the local atoms, so we don't need row and column
1276	>	! lookups for the types
1277	>
1278	>	me_i = atid(i)
1279	>
1280	>	! is the atom electrostatic?  See if it would have an
1281	>	! electrostatic interaction with itself
1282	>	iHash = InteractionHash(me_i,me_i)
1283
1284	<	if (FF_uses_RF .and. SIM_uses_RF) then
979	<
1284	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1285		#ifdef IS_MPI
1286	<	call scatter(rf_Row,rf,plan_atom_row_3d)
1287	<	call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
983	<	do i = 1,nlocal
984	<	rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
985	<	end do
986	<	#endif
987	<
988	<	do i = 1, nLocal
989	<
990	<	rfpot = 0.0_DP
991	<	#ifdef IS_MPI
992	<	me_i = atid_row(i)
1286	>	call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1287	>	t, do_pot)
1288		#else
1289	<	me_i = atid(i)
1289	>	call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1290	>	t, do_pot)
1291		#endif
1292	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1292	>	endif
1293	>
1294	>
1295	>	if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1296
1297	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1298	<
1299	<	mu_i = getDipoleMoment(me_i)
1300	<
1301	<	!! The reaction field needs to include a self contribution
1302	<	!! to the field:
1303	<	call accumulate_self_rf(i, mu_i, eFrame)
1304	<	!! Get the reaction field contribution to the
1305	<	!! potential and torques:
1306	<	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1297	>	! loop over the excludes to accumulate RF stuff we've
1298	>	! left out of the normal pair loop
1299	>
1300	>	do i1 = 1, nSkipsForAtom(i)
1301	>	j = skipsForAtom(i, i1)
1302	>
1303	>	! prevent overcounting of the skips
1304	>	if (i.lt.j) then
1305	>	call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1306	>	rVal = sqrt(ratmsq)
1307	>	call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1308		#ifdef IS_MPI
1309	<	pot_local = pot_local + rfpot
1309	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1310	>	vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1311		#else
1312	<	pot = pot + rfpot
1312	>	call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1313	>	vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1314	>	#endif
1315	>	endif
1316	>	enddo
1317	>	endif
1318
1319	+	if (do_box_dipole) then
1320	+	#ifdef IS_MPI
1321	+	call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1322	+	nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1323	+	pChgCount_local, nChgCount_local)
1324	+	#else
1325	+	call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1326	+	pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1327		#endif
1328	<	endif
1329	<	enddo
1016	<	endif
1328	>	endif
1329	>	enddo
1330		endif
1331
1019	–
1332		#ifdef IS_MPI
1021	–
1333		if (do_pot) then
1334	<	pot = pot + pot_local
1335	<	!! we assume the c code will do the allreduce to get the total potential
1336	<	!! we could do it right here if we needed to...
1334	>	#ifdef SINGLE_PRECISION
1335	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1336	>	mpi_comm_world,mpi_err)
1337	>	#else
1338	>	call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1339	>	mpi_sum, mpi_comm_world,mpi_err)
1340	>	#endif
1341		endif
1342	+
1343	+	if (do_box_dipole) then
1344
1345	<	if (do_stress) then
1346	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1347	<	mpi_comm_world,mpi_err)
1348	<	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1349	<	mpi_comm_world,mpi_err)
1350	<	endif
1351	<
1345	>	#ifdef SINGLE_PRECISION
1346	>	call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1347	>	mpi_comm_world, mpi_err)
1348	>	call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1349	>	mpi_comm_world, mpi_err)
1350	>	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1351	>	mpi_comm_world, mpi_err)
1352	>	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1353	>	mpi_comm_world, mpi_err)
1354	>	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1355	>	mpi_comm_world, mpi_err)
1356	>	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1357	>	mpi_comm_world, mpi_err)
1358	>	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1359	>	mpi_comm_world, mpi_err)
1360		#else
1361	+	call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1362	+	mpi_comm_world, mpi_err)
1363	+	call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1364	+	mpi_comm_world, mpi_err)
1365	+	call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1366	+	mpi_sum, mpi_comm_world, mpi_err)
1367	+	call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1368	+	mpi_sum, mpi_comm_world, mpi_err)
1369	+	call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1370	+	mpi_sum, mpi_comm_world, mpi_err)
1371	+	call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1372	+	mpi_sum, mpi_comm_world, mpi_err)
1373	+	call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1374	+	mpi_sum, mpi_comm_world, mpi_err)
1375	+	#endif
1376
1037	–	if (do_stress) then
1038	–	tau = tau_Temp
1039	–	virial = virial_Temp
1377		endif
1378	<
1378	>
1379		#endif
1380
1381	+	if (do_box_dipole) then
1382	+	! first load the accumulated dipole moment (if dipoles were present)
1383	+	boxDipole(1) = dipVec(1)
1384	+	boxDipole(2) = dipVec(2)
1385	+	boxDipole(3) = dipVec(3)
1386	+
1387	+	! now include the dipole moment due to charges
1388	+	! use the lesser of the positive and negative charge totals
1389	+	if (nChg .le. pChg) then
1390	+	chg_value = nChg
1391	+	else
1392	+	chg_value = pChg
1393	+	endif
1394	+
1395	+	! find the average positions
1396	+	if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1397	+	pChgPos = pChgPos / pChgCount
1398	+	nChgPos = nChgPos / nChgCount
1399	+	endif
1400	+
1401	+	! dipole is from the negative to the positive (physics notation)
1402	+	chgVec(1) = pChgPos(1) - nChgPos(1)
1403	+	chgVec(2) = pChgPos(2) - nChgPos(2)
1404	+	chgVec(3) = pChgPos(3) - nChgPos(3)
1405	+
1406	+	boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1407	+	boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1408	+	boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1409	+
1410	+	endif
1411	+
1412		end subroutine do_force_loop
1413
1414		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1415	<	eFrame, A, f, t, pot, vpair, fpair)
1415	>	eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1416
1417	<	real( kind = dp ) :: pot, vpair, sw
1417	>	real( kind = dp ) :: vpair, sw
1418	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1419		real( kind = dp ), dimension(3) :: fpair
1420		real( kind = dp ), dimension(nLocal)   :: mfact
1421		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1057 | Line 1426 | contains
1426		logical, intent(inout) :: do_pot
1427		integer, intent(in) :: i, j
1428		real ( kind = dp ), intent(inout) :: rijsq
1429	<	real ( kind = dp )                :: r
1429	>	real ( kind = dp ), intent(inout) :: r_grp
1430		real ( kind = dp ), intent(inout) :: d(3)
1431	<	real ( kind = dp ) :: ebalance
1431	>	real ( kind = dp ), intent(inout) :: d_grp(3)
1432	>	real ( kind = dp ), intent(inout) :: rCut
1433	>	real ( kind = dp ) :: r
1434	>	real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1435		integer :: me_i, me_j
1436	+	integer :: k
1437
1438	<	integer :: iMap
1438	>	integer :: iHash
1439
1440		r = sqrt(rijsq)
1441	<	vpair = 0.0d0
1442	<	fpair(1:3) = 0.0d0
1441	>
1442	>	vpair = 0.0_dp
1443	>	fpair(1:3) = 0.0_dp
1444
1445		#ifdef IS_MPI
1446		me_i = atid_row(i)
#	Line 1076 | Line 1450 | contains
1450		me_j = atid(j)
1451		#endif
1452
1453	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1454	<
1455	<	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1456	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1453	>	iHash = InteractionHash(me_i, me_j)
1454	>
1455	>	if ( iand(iHash, LJ_PAIR).ne.0 ) then
1456	>	call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1457	>	pot(VDW_POT), f, do_pot)
1458		endif
1459	<
1460	<	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1461	<	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1462	<	pot, eFrame, f, t, do_pot)
1088	<
1089	<	if (FF_uses_RF .and. SIM_uses_RF) then
1090	<
1091	<	! CHECK ME (RF needs to know about all electrostatic types)
1092	<	call accumulate_rf(i, j, r, eFrame, sw)
1093	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094	<	endif
1095	<
1459	>
1460	>	if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1461	>	call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1462	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1463		endif
1464	<
1465	<	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1464	>
1465	>	if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1466		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1467	<	pot, A, f, t, do_pot)
1467	>	pot(HB_POT), A, f, t, do_pot)
1468		endif
1469	<
1470	<	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1469	>
1470	>	if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1471		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1472	<	pot, A, f, t, do_pot)
1472	>	pot(HB_POT), A, f, t, do_pot)
1473		endif
1107	–
1108	–	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1109	–	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1110	–	pot, A, f, t, do_pot)
1111	–	endif
1474
1475	<	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1476	<	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1477	<	!           pot, A, f, t, do_pot)
1475	>	if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1476	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1477	>	pot(VDW_POT), A, f, t, do_pot)
1478		endif
1479	<
1480	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1481	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1482	<	do_pot)
1479	>
1480	>	if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1481	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1482	>	pot(VDW_POT), A, f, t, do_pot)
1483		endif
1484	<
1485	<	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1486	<	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1487	<	pot, A, f, t, do_pot)
1484	>
1485	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1486	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1487	>	pot(METALLIC_POT), f, do_pot)
1488		endif
1489	<
1490	<	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1489	>
1490	>	if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1491		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1492	<	pot, A, f, t, do_pot)
1492	>	pot(VDW_POT), A, f, t, do_pot)
1493		endif
1494
1495	+	if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1496	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1497	+	pot(VDW_POT), A, f, t, do_pot)
1498	+	endif
1499	+
1500	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1501	+	call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1502	+	pot(METALLIC_POT), f, do_pot)
1503	+	endif
1504	+
1505		end subroutine do_pair
1506
1507	<	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1507	>	subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1508		do_pot, do_stress, eFrame, A, f, t, pot)
1509
1510	<	real( kind = dp ) :: pot, sw
1510	>	real( kind = dp ) :: sw
1511	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1512		real( kind = dp ), dimension(9,nLocal) :: eFrame
1513		real (kind=dp), dimension(9,nLocal) :: A
1514		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1143 | Line 1516 | contains
1516
1517		logical, intent(inout) :: do_pot, do_stress
1518		integer, intent(in) :: i, j
1519	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1519	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1520		real ( kind = dp )                :: r, rc
1521		real ( kind = dp ), intent(inout) :: d(3), dc(3)
1522
1523	<	integer :: me_i, me_j, iMap
1523	>	integer :: me_i, me_j, iHash
1524
1525	+	r = sqrt(rijsq)
1526	+
1527		#ifdef IS_MPI
1528		me_i = atid_row(i)
1529		me_j = atid_col(j)
#	Line 1157 | Line 1532 | contains
1532		me_j = atid(j)
1533		#endif
1534
1535	<	iMap = InteractionMap(me_i, me_j)%InteractionHash
1535	>	iHash = InteractionHash(me_i, me_j)
1536
1537	<	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1538	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1537	>	if ( iand(iHash, EAM_PAIR).ne.0 ) then
1538	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1539		endif
1540	+
1541	+	if ( iand(iHash, SC_PAIR).ne.0 ) then
1542	+	call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1543	+	endif
1544
1545		end subroutine do_prepair
1546
1547
1548		subroutine do_preforce(nlocal,pot)
1549		integer :: nlocal
1550	<	real( kind = dp ) :: pot
1550	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1551
1552		if (FF_uses_EAM .and. SIM_uses_EAM) then
1553	<	call calc_EAM_preforce_Frho(nlocal,pot)
1553	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1554		endif
1555	<
1556	<
1555	>	if (FF_uses_SC .and. SIM_uses_SC) then
1556	>	call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1557	>	endif
1558		end subroutine do_preforce
1559
1560
#	Line 1186 | Line 1566 | contains
1566		real( kind = dp ) :: d(3), scaled(3)
1567		integer i
1568
1569	<	d(1:3) = q_j(1:3) - q_i(1:3)
1569	>	d(1) = q_j(1) - q_i(1)
1570	>	d(2) = q_j(2) - q_i(2)
1571	>	d(3) = q_j(3) - q_i(3)
1572
1573		! Wrap back into periodic box if necessary
1574		if ( SIM_uses_PBC ) then
1575
1576		if( .not.boxIsOrthorhombic ) then
1577		! calc the scaled coordinates.
1578	+	! scaled = matmul(HmatInv, d)
1579
1580	<	scaled = matmul(HmatInv, d)
1581	<
1580	>	scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1581	>	scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1582	>	scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1583	>
1584		! wrap the scaled coordinates
1585
1586	<	scaled = scaled  - anint(scaled)
1587	<
1586	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1587	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1588	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1589
1590		! calc the wrapped real coordinates from the wrapped scaled
1591		! coordinates
1592	+	! d = matmul(Hmat,scaled)
1593	+	d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1594	+	d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1595	+	d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1596
1207	–	d = matmul(Hmat,scaled)
1208	–
1597		else
1598		! calc the scaled coordinates.
1599
1600	<	do i = 1, 3
1601	<	scaled(i) = d(i) * HmatInv(i,i)
1600	>	scaled(1) = d(1) * HmatInv(1,1)
1601	>	scaled(2) = d(2) * HmatInv(2,2)
1602	>	scaled(3) = d(3) * HmatInv(3,3)
1603	>
1604	>	! wrap the scaled coordinates
1605	>
1606	>	scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1607	>	scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1608	>	scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1609
1610	<	! wrap the scaled coordinates
1610	>	! calc the wrapped real coordinates from the wrapped scaled
1611	>	! coordinates
1612
1613	<	scaled(i) = scaled(i) - anint(scaled(i))
1613	>	d(1) = scaled(1)*Hmat(1,1)
1614	>	d(2) = scaled(2)*Hmat(2,2)
1615	>	d(3) = scaled(3)*Hmat(3,3)
1616
1219	–	! calc the wrapped real coordinates from the wrapped scaled
1220	–	! coordinates
1221	–
1222	–	d(i) = scaled(i)*Hmat(i,i)
1223	–	enddo
1617		endif
1618
1619		endif
1620
1621	<	r_sq = dot_product(d,d)
1621	>	r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1622
1623		end subroutine get_interatomic_vector
1624
#	Line 1257 | Line 1650 | contains
1650		pot_Col = 0.0_dp
1651		pot_Temp = 0.0_dp
1652
1260	–	rf_Row = 0.0_dp
1261	–	rf_Col = 0.0_dp
1262	–	rf_Temp = 0.0_dp
1263	–
1653		#endif
1654
1655		if (FF_uses_EAM .and. SIM_uses_EAM) then
1656		call clean_EAM()
1657		endif
1658
1270	–	rf = 0.0_dp
1271	–	tau_Temp = 0.0_dp
1272	–	virial_Temp = 0.0_dp
1659		end subroutine zero_work_arrays
1660
1661		function skipThisPair(atom1, atom2) result(skip_it)
#	Line 1356 | Line 1742 | contains
1742
1743		function FF_UsesDirectionalAtoms() result(doesit)
1744		logical :: doesit
1745	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1361	<	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1745	>	doesit = FF_uses_DirectionalAtoms
1746		end function FF_UsesDirectionalAtoms
1747
1748		function FF_RequiresPrepairCalc() result(doesit)
1749		logical :: doesit
1750	<	doesit = FF_uses_EAM
1750	>	doesit = FF_uses_EAM .or. FF_uses_SC &
1751	>	.or. FF_uses_MEAM
1752		end function FF_RequiresPrepairCalc
1753
1369	–	function FF_RequiresPostpairCalc() result(doesit)
1370	–	logical :: doesit
1371	–	doesit = FF_uses_RF
1372	–	end function FF_RequiresPostpairCalc
1373	–
1754		#ifdef PROFILE
1755		function getforcetime() result(totalforcetime)
1756		real(kind=dp) :: totalforcetime
#	Line 1380 | Line 1760 | contains
1760
1761		!! This cleans componets of force arrays belonging only to fortran
1762
1763	<	subroutine add_stress_tensor(dpair, fpair)
1763	>	subroutine add_stress_tensor(dpair, fpair, tau)
1764
1765		real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1766	+	real( kind = dp ), dimension(9), intent(inout) :: tau
1767
1768		! because the d vector is the rj - ri vector, and
1769		! because fx, fy, fz are the force on atom i, we need a
1770		! negative sign here:
1771
1772	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1773	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1774	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1775	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1776	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1777	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1778	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1779	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1780	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1772	>	tau(1) = tau(1) - dpair(1) * fpair(1)
1773	>	tau(2) = tau(2) - dpair(1) * fpair(2)
1774	>	tau(3) = tau(3) - dpair(1) * fpair(3)
1775	>	tau(4) = tau(4) - dpair(2) * fpair(1)
1776	>	tau(5) = tau(5) - dpair(2) * fpair(2)
1777	>	tau(6) = tau(6) - dpair(2) * fpair(3)
1778	>	tau(7) = tau(7) - dpair(3) * fpair(1)
1779	>	tau(8) = tau(8) - dpair(3) * fpair(2)
1780	>	tau(9) = tau(9) - dpair(3) * fpair(3)
1781
1401	–	virial_Temp = virial_Temp + &
1402	–	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1403	–
1782		end subroutine add_stress_tensor
1783
1784		end module doForces

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines