ViewVC Help

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

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

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 1610 by gezelter, Wed Oct 20 04:19:55 2004 UTC vs.
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC

#	Line 1 | Line 1
1	+	!!
2	+	!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	+	!!
4	+	!! The University of Notre Dame grants you ("Licensee") a
5	+	!! non-exclusive, royalty free, license to use, modify and
6	+	!! redistribute this software in source and binary code form, provided
7	+	!! that the following conditions are met:
8	+	!!
9	+	!! 1. Acknowledgement of the program authors must be made in any
10	+	!!    publication of scientific results based in part on use of the
11	+	!!    program.  An acceptable form of acknowledgement is citation of
12	+	!!    the article in which the program was described (Matthew
13	+	!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	+	!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	+	!!    Parallel Simulation Engine for Molecular Dynamics,"
16	+	!!    J. Comput. Chem. 26, pp. 252-271 (2005))
17	+	!!
18	+	!! 2. Redistributions of source code must retain the above copyright
19	+	!!    notice, this list of conditions and the following disclaimer.
20	+	!!
21	+	!! 3. Redistributions in binary form must reproduce the above copyright
22	+	!!    notice, this list of conditions and the following disclaimer in the
23	+	!!    documentation and/or other materials provided with the
24	+	!!    distribution.
25	+	!!
26	+	!! This software is provided "AS IS," without a warranty of any
27	+	!! kind. All express or implied conditions, representations and
28	+	!! warranties, including any implied warranty of merchantability,
29	+	!! fitness for a particular purpose or non-infringement, are hereby
30	+	!! excluded.  The University of Notre Dame and its licensors shall not
31	+	!! be liable for any damages suffered by licensee as a result of
32	+	!! using, modifying or distributing the software or its
33	+	!! derivatives. In no event will the University of Notre Dame or its
34	+	!! licensors be liable for any lost revenue, profit or data, or for
35	+	!! direct, indirect, special, consequential, incidental or punitive
36	+	!! damages, however caused and regardless of the theory of liability,
37	+	!! arising out of the use of or inability to use software, even if the
38	+	!! University of Notre Dame has been advised of the possibility of
39	+	!! such damages.
40	+	!!
41	+
42		!! doForces.F90
43		!! module doForces
44		!! Calculates Long Range forces.
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.1 2004-10-20 04:19:55 gezelter Exp $, $Date: 2004-10-20 04:19:55 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
48	>	!! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
49
50	+
51		module doForces
52		use force_globals
53		use simulation
#	Line 14 | Line 56 | module doForces
56		use switcheroo
57		use neighborLists
58		use lj
59	<	use sticky_pair
60	<	use dipole_dipole
19	<	use charge_charge
59	>	use sticky
60	>	use electrostatic_module
61		use reaction_field
62		use gb_pair
63	+	use shapes
64		use vector_class
65		use eam
66		use status
#	Line 30 | Line 72 | module doForces
72		PRIVATE
73
74		#define __FORTRAN90
33	–	#include "UseTheForce/fForceField.h"
75		#include "UseTheForce/fSwitchingFunction.h"
76	+	#include "UseTheForce/DarkSide/fInteractionMap.h"
77
78		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
79		INTEGER, PARAMETER:: PAIR_LOOP    = 2
80
81		logical, save :: haveRlist = .false.
82		logical, save :: haveNeighborList = .false.
41	–	logical, save :: havePolicies = .false.
83		logical, save :: haveSIMvariables = .false.
43	–	logical, save :: havePropertyMap = .false.
84		logical, save :: haveSaneForceField = .false.
85	<	logical, save :: FF_uses_LJ
86	<	logical, save :: FF_uses_sticky
87	<	logical, save :: FF_uses_charges
88	<	logical, save :: FF_uses_dipoles
89	<	logical, save :: FF_uses_RF
90	<	logical, save :: FF_uses_GB
85	>	logical, save :: haveInteractionMap = .false.
86	>
87	>	logical, save :: FF_uses_DirectionalAtoms
88	>	logical, save :: FF_uses_LennardJones
89	>	logical, save :: FF_uses_Electrostatics
90	>	logical, save :: FF_uses_Charges
91	>	logical, save :: FF_uses_Dipoles
92	>	logical, save :: FF_uses_Quadrupoles
93	>	logical, save :: FF_uses_Sticky
94	>	logical, save :: FF_uses_StickyPower
95	>	logical, save :: FF_uses_GayBerne
96		logical, save :: FF_uses_EAM
97	<	logical, save :: SIM_uses_LJ
98	<	logical, save :: SIM_uses_sticky
99	<	logical, save :: SIM_uses_charges
100	<	logical, save :: SIM_uses_dipoles
101	<	logical, save :: SIM_uses_RF
102	<	logical, save :: SIM_uses_GB
97	>	logical, save :: FF_uses_Shapes
98	>	logical, save :: FF_uses_FLARB
99	>	logical, save :: FF_uses_RF
100	>
101	>	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
110		logical, save :: SIM_uses_EAM
111	+	logical, save :: SIM_uses_Shapes
112	+	logical, save :: SIM_uses_FLARB
113	+	logical, save :: SIM_uses_RF
114		logical, save :: SIM_requires_postpair_calc
115		logical, save :: SIM_requires_prepair_calc
61	–	logical, save :: SIM_uses_directional_atoms
116		logical, save :: SIM_uses_PBC
117		logical, save :: SIM_uses_molecular_cutoffs
118
119	<	real(kind=dp), save :: rlist, rlistsq
119	>	!!!GO AWAY---------
120	>	!!!!!real(kind=dp), save :: rlist, rlistsq
121
122		public :: init_FF
123		public :: do_force_loop
124	<	public :: setRlistDF
124	>	!  public :: setRlistDF
125	>	!public :: addInteraction
126	>	!public :: setInteractionHash
127	>	!public :: getInteractionHash
128	>	public :: createInteractionMap
129	>	public :: createRcuts
130
131		#ifdef PROFILE
132		public :: getforcetime
#	Line 75 | Line 135 | module doForces
135		integer :: nLoops
136		#endif
137
138	<	type :: Properties
139	<	logical :: is_lj     = .false.
140	<	logical :: is_sticky = .false.
141	<	logical :: is_dp     = .false.
142	<	logical :: is_gb     = .false.
143	<	logical :: is_eam    = .false.
144	<	logical :: is_charge = .false.
145	<	real(kind=DP) :: charge = 0.0_DP
86	<	real(kind=DP) :: dipole_moment = 0.0_DP
87	<	end type Properties
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
143	>
144	>	type(Interaction), dimension(:,:),allocatable :: InteractionMap
145	>
146
147	<	type(Properties), dimension(:),allocatable :: PropertyMap
90	<
147	>
148		contains
149
150	<	subroutine setRlistDF( this_rlist )
151	<
95	<	real(kind=dp) :: this_rlist
96	<
97	<	rlist = this_rlist
98	<	rlistsq = rlist * rlist
99	<
100	<	haveRlist = .true.
101	<
102	<	end subroutine setRlistDF
103	<
104	<	subroutine createPropertyMap(status)
150	>
151	>	subroutine createInteractionMap(status)
152		integer :: nAtypes
153	<	integer :: status
153	>	integer, intent(out) :: status
154		integer :: i
155	<	logical :: thisProperty
156	<	real (kind=DP) :: thisDPproperty
157	<
155	>	integer :: j
156	>	integer :: ihash
157	>	real(kind=dp) :: myRcut
158	>	!! Test Types
159	>	logical :: i_is_LJ
160	>	logical :: i_is_Elect
161	>	logical :: i_is_Sticky
162	>	logical :: i_is_StickyP
163	>	logical :: i_is_GB
164	>	logical :: i_is_EAM
165	>	logical :: i_is_Shape
166	>	logical :: j_is_LJ
167	>	logical :: j_is_Elect
168	>	logical :: j_is_Sticky
169	>	logical :: j_is_StickyP
170	>	logical :: j_is_GB
171	>	logical :: j_is_EAM
172	>	logical :: j_is_Shape
173	>
174		status = 0
175	<
175	>
176	>	if (.not. associated(atypes)) then
177	>	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
178	>	status = -1
179	>	return
180	>	endif
181	>
182		nAtypes = getSize(atypes)
183	<
183	>
184		if (nAtypes == 0) then
185		status = -1
186		return
187		end if
119	–
120	–	if (.not. allocated(PropertyMap)) then
121	–	allocate(PropertyMap(nAtypes))
122	–	endif
188
189	+	if (.not. allocated(InteractionMap)) then
190	+	allocate(InteractionMap(nAtypes,nAtypes))
191	+	endif
192	+
193		do i = 1, nAtypes
194	<	call getElementProperty(atypes, i, "is_LJ", thisProperty)
195	<	PropertyMap(i)%is_LJ = thisProperty
194	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
195	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
196	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
197	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
198	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
199	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
200	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
201
202	<	call getElementProperty(atypes, i, "is_Charge", thisProperty)
129	<	PropertyMap(i)%is_Charge = thisProperty
130	<
131	<	if (thisProperty) then
132	<	call getElementProperty(atypes, i, "charge", thisDPproperty)
133	<	PropertyMap(i)%charge = thisDPproperty
134	<	endif
202	>	do j = i, nAtypes
203
204	<	call getElementProperty(atypes, i, "is_DP", thisProperty)
205	<	PropertyMap(i)%is_DP = thisProperty
204	>	iHash = 0
205	>	myRcut = 0.0_dp
206
207	<	if (thisProperty) then
208	<	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
209	<	PropertyMap(i)%dipole_moment = thisDPproperty
210	<	endif
207	>	call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
208	>	call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
209	>	call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
210	>	call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
211	>	call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
212	>	call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
213	>	call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
214
215	<	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
216	<	PropertyMap(i)%is_Sticky = thisProperty
217	<	call getElementProperty(atypes, i, "is_GB", thisProperty)
218	<	PropertyMap(i)%is_GB = thisProperty
219	<	call getElementProperty(atypes, i, "is_EAM", thisProperty)
220	<	PropertyMap(i)%is_EAM = thisProperty
215	>	if (i_is_LJ .and. j_is_LJ) then
216	>	iHash = ior(iHash, LJ_PAIR)
217	>	endif
218	>
219	>	if (i_is_Elect .and. j_is_Elect) then
220	>	iHash = ior(iHash, ELECTROSTATIC_PAIR)
221	>	endif
222	>
223	>	if (i_is_Sticky .and. j_is_Sticky) then
224	>	iHash = ior(iHash, STICKY_PAIR)
225	>	endif
226	>
227	>	if (i_is_StickyP .and. j_is_StickyP) then
228	>	iHash = ior(iHash, STICKYPOWER_PAIR)
229	>	endif
230	>
231	>	if (i_is_EAM .and. j_is_EAM) then
232	>	iHash = ior(iHash, EAM_PAIR)
233	>	endif
234	>
235	>	if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
236	>	if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
237	>	if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
238	>
239	>	if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
240	>	if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
241	>	if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
242	>
243	>
244	>	InteractionMap(i,j)%InteractionHash = iHash
245	>	InteractionMap(j,i)%InteractionHash = iHash
246	>
247	>	end do
248	>
249		end do
250
251	<	havePropertyMap = .true.
251	>	haveInteractionMap = .true.
252	>	end subroutine createInteractionMap
253
254	<	end subroutine createPropertyMap
254	>	! 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.
255	>	subroutine createRcuts(defaultRList,stat)
256	>	real(kind=dp), intent(in), optional :: defaultRList
257	>	integer :: iMap
258	>	integer :: map_i,map_j
259	>	real(kind=dp) :: thisRCut = 0.0_dp
260	>	real(kind=dp) :: actualCutoff = 0.0_dp
261	>	integer, intent(out) :: stat
262	>	integer :: nAtypes
263	>	integer :: myStatus
264
265	+	stat = 0
266	+	if (.not. haveInteractionMap) then
267	+
268	+	call createInteractionMap(myStatus)
269	+
270	+	if (myStatus .ne. 0) then
271	+	write(default_error, *) 'createInteractionMap failed in doForces!'
272	+	stat = -1
273	+	return
274	+	endif
275	+	endif
276	+
277	+
278	+	nAtypes = getSize(atypes)
279	+	!! If we pass a default rcut, set all atypes to that cutoff distance
280	+	if(present(defaultRList)) then
281	+	InteractionMap(:,:)%rList = defaultRList
282	+	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
283	+	haveRlist = .true.
284	+	return
285	+	end if
286	+
287	+	do map_i = 1,nAtypes
288	+	do map_j = map_i,nAtypes
289	+	iMap = InteractionMap(map_i, map_j)%InteractionHash
290	+
291	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
292	+	! thisRCut = getLJCutOff(map_i,map_j)
293	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
294	+	endif
295	+
296	+	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
297	+	! thisRCut = getElectrostaticCutOff(map_i,map_j)
298	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
299	+	endif
300	+
301	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
302	+	! thisRCut = getStickyCutOff(map_i,map_j)
303	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
304	+	endif
305	+
306	+	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307	+	! thisRCut = getStickyPowerCutOff(map_i,map_j)
308	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
309	+	endif
310	+
311	+	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312	+	! thisRCut = getGayberneCutOff(map_i,map_j)
313	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
314	+	endif
315	+
316	+	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317	+	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
319	+	endif
320	+
321	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
322	+	!              thisRCut = getEAMCutOff(map_i,map_j)
323	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
324	+	endif
325	+
326	+	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
327	+	!              thisRCut = getShapeCutOff(map_i,map_j)
328	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
329	+	endif
330	+
331	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
332	+	!              thisRCut = getShapeLJCutOff(map_i,map_j)
333	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
334	+	endif
335	+	InteractionMap(map_i, map_j)%rList = actualCutoff
336	+	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337	+	end do
338	+	end do
339	+	haveRlist = .true.
340	+	end subroutine createRcuts
341	+
342	+
343	+	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
344	+	!!$  subroutine setRlistDF( this_rlist )
345	+	!!$
346	+	!!$   real(kind=dp) :: this_rlist
347	+	!!$
348	+	!!$    rlist = this_rlist
349	+	!!$    rlistsq = rlist * rlist
350	+	!!$
351	+	!!$    haveRlist = .true.
352	+	!!$
353	+	!!$  end subroutine setRlistDF
354	+
355	+
356		subroutine setSimVariables()
357	<	SIM_uses_LJ = SimUsesLJ()
358	<	SIM_uses_sticky = SimUsesSticky()
359	<	SIM_uses_charges = SimUsesCharges()
360	<	SIM_uses_dipoles = SimUsesDipoles()
361	<	SIM_uses_RF = SimUsesRF()
362	<	SIM_uses_GB = SimUsesGB()
357	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358	>	SIM_uses_LennardJones = SimUsesLennardJones()
359	>	SIM_uses_Electrostatics = SimUsesElectrostatics()
360	>	SIM_uses_Charges = SimUsesCharges()
361	>	SIM_uses_Dipoles = SimUsesDipoles()
362	>	SIM_uses_Sticky = SimUsesSticky()
363	>	SIM_uses_StickyPower = SimUsesStickyPower()
364	>	SIM_uses_GayBerne = SimUsesGayBerne()
365		SIM_uses_EAM = SimUsesEAM()
366	+	SIM_uses_Shapes = SimUsesShapes()
367	+	SIM_uses_FLARB = SimUsesFLARB()
368	+	SIM_uses_RF = SimUsesRF()
369		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
370		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
166	–	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
371		SIM_uses_PBC = SimUsesPBC()
168	–	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
372
373		haveSIMvariables = .true.
374
#	Line 178 | Line 381 | contains
381		integer :: myStatus
382
383		error = 0
384	<
385	<	if (.not. havePropertyMap) then
384	>
385	>	if (.not. haveInteractionMap) then
386
387		myStatus = 0
388
389	<	call createPropertyMap(myStatus)
389	>	call createInteractionMap(myStatus)
390
391		if (myStatus .ne. 0) then
392	<	write(default_error, *) 'createPropertyMap failed in doForces!'
392	>	write(default_error, *) 'createInteractionMap failed in doForces!'
393		error = -1
394		return
395		endif
#	Line 202 | Line 405 | contains
405		return
406		endif
407
205	–	if (SIM_uses_LJ .and. FF_uses_LJ) then
206	–	if (.not. havePolicies) then
207	–	write(default_error, *) 'LJ mixing Policies have not been set in doForces!'
208	–	error = -1
209	–	return
210	–	endif
211	–	endif
212	–
408		if (.not. haveNeighborList) then
409		write(default_error, *) 'neighbor list has not been initialized in doForces!'
410		error = -1
#	Line 231 | Line 426 | contains
426		#endif
427		return
428		end subroutine doReadyCheck
234	–
429
236	–	subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
430
431	<	integer, intent(in) :: LJMIXPOLICY
431	>	subroutine init_FF(use_RF_c, thisStat)
432	>
433		logical, intent(in) :: use_RF_c
434
435		integer, intent(out) :: thisStat
#	Line 248 | Line 442 | contains
442
443		!! Fortran's version of a cast:
444		FF_uses_RF = use_RF_c
445	<
445	>
446		!! init_FF is called *after* all of the atom types have been
447		!! defined in atype_module using the new_atype subroutine.
448		!!
449		!! this will scan through the known atypes and figure out what
450		!! interactions are used by the force field.
451	<
452	<	FF_uses_LJ = .false.
453	<	FF_uses_sticky = .false.
454	<	FF_uses_charges = .false.
455	<	FF_uses_dipoles = .false.
456	<	FF_uses_GB = .false.
451	>
452	>	FF_uses_DirectionalAtoms = .false.
453	>	FF_uses_LennardJones = .false.
454	>	FF_uses_Electrostatics = .false.
455	>	FF_uses_Charges = .false.
456	>	FF_uses_Dipoles = .false.
457	>	FF_uses_Sticky = .false.
458	>	FF_uses_StickyPower = .false.
459	>	FF_uses_GayBerne = .false.
460		FF_uses_EAM = .false.
461	<
462	<	call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
266	<	if (nMatches .gt. 0) FF_uses_LJ = .true.
461	>	FF_uses_Shapes = .false.
462	>	FF_uses_FLARB = .false.
463
464	<	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
465	<	if (nMatches .gt. 0) FF_uses_charges = .true.
464	>	call getMatchingElementList(atypes, "is_Directional", .true., &
465	>	nMatches, MatchList)
466	>	if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
467
468	<	call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
469	<	if (nMatches .gt. 0) FF_uses_dipoles = .true.
470	<
468	>	call getMatchingElementList(atypes, "is_LennardJones", .true., &
469	>	nMatches, MatchList)
470	>	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471	>
472	>	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473	>	nMatches, MatchList)
474	>	if (nMatches .gt. 0) then
475	>	FF_uses_Electrostatics = .true.
476	>	endif
477	>
478	>	call getMatchingElementList(atypes, "is_Charge", .true., &
479	>	nMatches, MatchList)
480	>	if (nMatches .gt. 0) then
481	>	FF_uses_Charges = .true.
482	>	FF_uses_Electrostatics = .true.
483	>	endif
484	>
485	>	call getMatchingElementList(atypes, "is_Dipole", .true., &
486	>	nMatches, MatchList)
487	>	if (nMatches .gt. 0) then
488	>	FF_uses_Dipoles = .true.
489	>	FF_uses_Electrostatics = .true.
490	>	FF_uses_DirectionalAtoms = .true.
491	>	endif
492	>
493	>	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494	>	nMatches, MatchList)
495	>	if (nMatches .gt. 0) then
496	>	FF_uses_Quadrupoles = .true.
497	>	FF_uses_Electrostatics = .true.
498	>	FF_uses_DirectionalAtoms = .true.
499	>	endif
500	>
501		call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502		MatchList)
503	<	if (nMatches .gt. 0) FF_uses_Sticky = .true.
504	<
505	<	call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
506	<	if (nMatches .gt. 0) FF_uses_GB = .true.
503	>	if (nMatches .gt. 0) then
504	>	FF_uses_Sticky = .true.
505	>	FF_uses_DirectionalAtoms = .true.
506	>	endif
507	>
508	>	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509	>	MatchList)
510	>	if (nMatches .gt. 0) then
511	>	FF_uses_StickyPower = .true.
512	>	FF_uses_DirectionalAtoms = .true.
513	>	endif
514
515	+	call getMatchingElementList(atypes, "is_GayBerne", .true., &
516	+	nMatches, MatchList)
517	+	if (nMatches .gt. 0) then
518	+	FF_uses_GayBerne = .true.
519	+	FF_uses_DirectionalAtoms = .true.
520	+	endif
521	+
522		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
523		if (nMatches .gt. 0) FF_uses_EAM = .true.
524	<
524	>
525	>	call getMatchingElementList(atypes, "is_Shape", .true., &
526	>	nMatches, MatchList)
527	>	if (nMatches .gt. 0) then
528	>	FF_uses_Shapes = .true.
529	>	FF_uses_DirectionalAtoms = .true.
530	>	endif
531	>
532	>	call getMatchingElementList(atypes, "is_FLARB", .true., &
533	>	nMatches, MatchList)
534	>	if (nMatches .gt. 0) FF_uses_FLARB = .true.
535	>
536		!! Assume sanity (for the sake of argument)
537		haveSaneForceField = .true.
538
539		!! check to make sure the FF_uses_RF setting makes sense
540	<
540	>
541		if (FF_uses_dipoles) then
542		if (FF_uses_RF) then
543		dielect = getDielect()
#	Line 298 | Line 550 | contains
550		haveSaneForceField = .false.
551		return
552		endif
301	–	endif
302	–
303	–	if (FF_uses_LJ) then
304	–
305	–	select case (LJMIXPOLICY)
306	–	case (LB_MIXING_RULE)
307	–	call init_lj_FF(LB_MIXING_RULE, my_status)
308	–	case (EXPLICIT_MIXING_RULE)
309	–	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
310	–	case default
311	–	write(default_error,*) 'unknown LJ Mixing Policy!'
312	–	thisStat = -1
313	–	haveSaneForceField = .false.
314	–	return
315	–	end select
316	–	if (my_status /= 0) then
317	–	thisStat = -1
318	–	haveSaneForceField = .false.
319	–	return
320	–	end if
321	–	havePolicies = .true.
553		endif
554
555	<	if (FF_uses_sticky) then
556	<	call check_sticky_FF(my_status)
557	<	if (my_status /= 0) then
558	<	thisStat = -1
559	<	haveSaneForceField = .false.
560	<	return
561	<	end if
562	<	endif
555	>	!sticky module does not contain check_sticky_FF anymore
556	>	!if (FF_uses_sticky) then
557	>	!   call check_sticky_FF(my_status)
558	>	!   if (my_status /= 0) then
559	>	!      thisStat = -1
560	>	!      haveSaneForceField = .false.
561	>	!      return
562	>	!   end if
563	>	!endif
564
333	–
565		if (FF_uses_EAM) then
566	<	call init_EAM_FF(my_status)
566	>	call init_EAM_FF(my_status)
567		if (my_status /= 0) then
568		write(default_error, *) "init_EAM_FF returned a bad status"
569		thisStat = -1
#	Line 341 | Line 572 | contains
572		end if
573		endif
574
575	<	if (FF_uses_GB) then
575	>	if (FF_uses_GayBerne) then
576		call check_gb_pair_FF(my_status)
577		if (my_status .ne. 0) then
578		thisStat = -1
#	Line 350 | Line 581 | contains
581		endif
582		endif
583
584	<	if (FF_uses_GB .and. FF_uses_LJ) then
584	>	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585		endif
586	+
587		if (.not. haveNeighborList) then
588		!! Create neighbor lists
589		call expandNeighborList(nLocal, my_status)
#	Line 363 | Line 595 | contains
595		haveNeighborList = .true.
596		endif
597
366	–
367	–
598		end subroutine init_FF
369	–
599
600	+
601		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
602		!------------------------------------------------------------->
603	<	subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
603	>	subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
604		do_pot_c, do_stress_c, error)
605		!! Position array provided by C, dimensioned by getNlocal
606		real ( kind = dp ), dimension(3, nLocal) :: q
#	Line 379 | Line 609 | contains
609		!! Rotation Matrix for each long range particle in simulation.
610		real( kind = dp), dimension(9, nLocal) :: A
611		!! Unit vectors for dipoles (lab frame)
612	<	real( kind = dp ), dimension(3,nLocal) :: u_l
612	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
613		!! Force array provided by C, dimensioned by getNlocal
614		real ( kind = dp ), dimension(3,nLocal) :: f
615		!! Torsion array provided by C, dimensioned by getNlocal
#	Line 417 | Line 647 | contains
647		integer :: localError
648		integer :: propPack_i, propPack_j
649		integer :: loopStart, loopEnd, loop
650	<
650	>	integer :: iMap
651		real(kind=dp) :: listSkin = 1.0
652	<
652	>
653		!! initialize local variables
654	<
654	>
655		#ifdef IS_MPI
656		pot_local = 0.0_dp
657		nAtomsInRow   = getNatomsInRow(plan_atom_row)
#	Line 431 | Line 661 | contains
661		#else
662		natoms = nlocal
663		#endif
664	<
664	>
665		call doReadyCheck(localError)
666		if ( localError .ne. 0 ) then
667		call handleError("do_force_loop", "Not Initialized")
#	Line 439 | Line 669 | contains
669		return
670		end if
671		call zero_work_arrays()
672	<
672	>
673		do_pot = do_pot_c
674		do_stress = do_stress_c
675	<
675	>
676		! Gather all information needed by all force loops:
677	<
677	>
678		#ifdef IS_MPI
679	<
679	>
680		call gather(q, q_Row, plan_atom_row_3d)
681		call gather(q, q_Col, plan_atom_col_3d)
682
683		call gather(q_group, q_group_Row, plan_group_row_3d)
684		call gather(q_group, q_group_Col, plan_group_col_3d)
685	<
686	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
687	<	call gather(u_l, u_l_Row, plan_atom_row_3d)
688	<	call gather(u_l, u_l_Col, plan_atom_col_3d)
689	<
685	>
686	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
687	>	call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
688	>	call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
689	>
690		call gather(A, A_Row, plan_atom_row_rotation)
691		call gather(A, A_Col, plan_atom_col_rotation)
692		endif
693	<
693	>
694		#endif
695	<
695	>
696		!! Begin force loop timing:
697		#ifdef PROFILE
698		call cpu_time(forceTimeInitial)
699		nloops = nloops + 1
700		#endif
701	<
701	>
702		loopEnd = PAIR_LOOP
703		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
704		loopStart = PREPAIR_LOOP
#	Line 483 | Line 713 | contains
713		if (loop .eq. loopStart) then
714		#ifdef IS_MPI
715		call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
716	<	update_nlist)
716	>	update_nlist)
717		#else
718		call checkNeighborList(nGroups, q_group, listSkin, &
719	<	update_nlist)
719	>	update_nlist)
720		#endif
721		endif
722	<
722	>
723		if (update_nlist) then
724		!! save current configuration and construct neighbor list
725		#ifdef IS_MPI
#	Line 500 | Line 730 | contains
730		neighborListSize = size(list)
731		nlist = 0
732		endif
733	<
733	>
734		istart = 1
735		#ifdef IS_MPI
736		iend = nGroupsInRow
#	Line 509 | Line 739 | contains
739		#endif
740		outer: do i = istart, iend
741
742	+	#ifdef IS_MPI
743	+	me_i = atid_row(i)
744	+	#else
745	+	me_i = atid(i)
746	+	#endif
747	+
748		if (update_nlist) point(i) = nlist + 1
749	<
749	>
750		n_in_i = groupStartRow(i+1) - groupStartRow(i)
751	<
751	>
752		if (update_nlist) then
753		#ifdef IS_MPI
754		jstart = 1
#	Line 527 | Line 763 | contains
763		! make sure group i has neighbors
764		if (jstart .gt. jend) cycle outer
765		endif
766	<
766	>
767		do jnab = jstart, jend
768		if (update_nlist) then
769		j = jnab
#	Line 536 | Line 772 | contains
772		endif
773
774		#ifdef IS_MPI
775	+	me_j = atid_col(j)
776		call get_interatomic_vector(q_group_Row(:,i), &
777		q_group_Col(:,j), d_grp, rgrpsq)
778		#else
779	+	me_j = atid(j)
780		call get_interatomic_vector(q_group(:,i), &
781		q_group(:,j), d_grp, rgrpsq)
782		#endif
783
784	<	if (rgrpsq < rlistsq) then
784	>	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
785		if (update_nlist) then
786		nlist = nlist + 1
787	<
787	>
788		if (nlist > neighborListSize) then
789		#ifdef IS_MPI
790		call expandNeighborList(nGroupsInRow, listerror)
#	Line 560 | Line 798 | contains
798		end if
799		neighborListSize = size(list)
800		endif
801	<
801	>
802		list(nlist) = j
803		endif
804	<
804	>
805		if (loop .eq. PAIR_LOOP) then
806		vij = 0.0d0
807		fij(1:3) = 0.0d0
808		endif
809	<
809	>
810		call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
811		in_switching_region)
812	<
812	>
813		n_in_j = groupStartCol(j+1) - groupStartCol(j)
814	<
814	>
815		do ia = groupStartRow(i), groupStartRow(i+1)-1
816	<
816	>
817		atom1 = groupListRow(ia)
818	<
818	>
819		inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
820	<
820	>
821		atom2 = groupListCol(jb)
822	<
822	>
823		if (skipThisPair(atom1, atom2)) cycle inner
824
825		if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
#	Line 601 | Line 839 | contains
839		#ifdef IS_MPI
840		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
841		rgrpsq, d_grp, do_pot, do_stress, &
842	<	u_l, A, f, t, pot_local)
842	>	eFrame, A, f, t, pot_local)
843		#else
844		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
845		rgrpsq, d_grp, do_pot, do_stress, &
846	<	u_l, A, f, t, pot)
846	>	eFrame, A, f, t, pot)
847		#endif
848		else
849		#ifdef IS_MPI
850		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
851		do_pot, &
852	<	u_l, A, f, t, pot_local, vpair, fpair)
852	>	eFrame, A, f, t, pot_local, vpair, fpair)
853		#else
854		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
855		do_pot,  &
856	<	u_l, A, f, t, pot, vpair, fpair)
856	>	eFrame, A, f, t, pot, vpair, fpair)
857		#endif
858
859		vij = vij + vpair
#	Line 623 | Line 861 | contains
861		endif
862		enddo inner
863		enddo
864	<
864	>
865		if (loop .eq. PAIR_LOOP) then
866		if (in_switching_region) then
867		swderiv = vij*dswdr/rgrp
868		fij(1) = fij(1) + swderiv*d_grp(1)
869		fij(2) = fij(2) + swderiv*d_grp(2)
870		fij(3) = fij(3) + swderiv*d_grp(3)
871	<
871	>
872		do ia=groupStartRow(i), groupStartRow(i+1)-1
873		atom1=groupListRow(ia)
874		mf = mfactRow(atom1)
#	Line 644 | Line 882 | contains
882		f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
883		#endif
884		enddo
885	<
885	>
886		do jb=groupStartCol(j), groupStartCol(j+1)-1
887		atom2=groupListCol(jb)
888		mf = mfactCol(atom2)
#	Line 659 | Line 897 | contains
897		#endif
898		enddo
899		endif
900	<
900	>
901		if (do_stress) call add_stress_tensor(d_grp, fij)
902		endif
903		end if
904		enddo
905		enddo outer
906	<
906	>
907		if (update_nlist) then
908		#ifdef IS_MPI
909		point(nGroupsInRow + 1) = nlist + 1
#	Line 679 | Line 917 | contains
917		update_nlist = .false.
918		endif
919		endif
920	<
920	>
921		if (loop .eq. PREPAIR_LOOP) then
922		call do_preforce(nlocal, pot)
923		endif
924	<
924	>
925		enddo
926	<
926	>
927		!! Do timing
928		#ifdef PROFILE
929		call cpu_time(forceTimeFinal)
930		forceTime = forceTime + forceTimeFinal - forceTimeInitial
931		#endif
932	<
932	>
933		#ifdef IS_MPI
934		!!distribute forces
935	<
935	>
936		f_temp = 0.0_dp
937		call scatter(f_Row,f_temp,plan_atom_row_3d)
938		do i = 1,nlocal
939		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
940		end do
941	<
941	>
942		f_temp = 0.0_dp
943		call scatter(f_Col,f_temp,plan_atom_col_3d)
944		do i = 1,nlocal
945		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
946		end do
947	<
948	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
947	>
948	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
949		t_temp = 0.0_dp
950		call scatter(t_Row,t_temp,plan_atom_row_3d)
951		do i = 1,nlocal
#	Line 715 | Line 953 | contains
953		end do
954		t_temp = 0.0_dp
955		call scatter(t_Col,t_temp,plan_atom_col_3d)
956	<
956	>
957		do i = 1,nlocal
958		t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
959		end do
960		endif
961	<
961	>
962		if (do_pot) then
963		! scatter/gather pot_row into the members of my column
964		call scatter(pot_Row, pot_Temp, plan_atom_row)
965	<
965	>
966		! scatter/gather pot_local into all other procs
967		! add resultant to get total pot
968		do i = 1, nlocal
969		pot_local = pot_local + pot_Temp(i)
970		enddo
971	<
971	>
972		pot_Temp = 0.0_DP
973	<
973	>
974		call scatter(pot_Col, pot_Temp, plan_atom_col)
975		do i = 1, nlocal
976		pot_local = pot_local + pot_Temp(i)
977		enddo
978	<
978	>
979		endif
980		#endif
981	<
981	>
982		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
983	<
983	>
984		if (FF_uses_RF .and. SIM_uses_RF) then
985	<
985	>
986		#ifdef IS_MPI
987		call scatter(rf_Row,rf,plan_atom_row_3d)
988		call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
#	Line 752 | Line 990 | contains
990		rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
991		end do
992		#endif
993	<
993	>
994		do i = 1, nLocal
995	<
995	>
996		rfpot = 0.0_DP
997		#ifdef IS_MPI
998		me_i = atid_row(i)
999		#else
1000		me_i = atid(i)
1001		#endif
1002	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
1003
1004	<	if (PropertyMap(me_i)%is_DP) then
1005	<
1006	<	mu_i = PropertyMap(me_i)%dipole_moment
1007	<
1004	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1005	>
1006	>	mu_i = getDipoleMoment(me_i)
1007	>
1008		!! The reaction field needs to include a self contribution
1009		!! to the field:
1010	<	call accumulate_self_rf(i, mu_i, u_l)
1010	>	call accumulate_self_rf(i, mu_i, eFrame)
1011		!! Get the reaction field contribution to the
1012		!! potential and torques:
1013	<	call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
1013	>	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1014		#ifdef IS_MPI
1015		pot_local = pot_local + rfpot
1016		#else
1017		pot = pot + rfpot
1018	<
1018	>
1019		#endif
1020	<	endif
1020	>	endif
1021		enddo
1022		endif
1023		endif
1024	<
1025	<
1024	>
1025	>
1026		#ifdef IS_MPI
1027	<
1027	>
1028		if (do_pot) then
1029		pot = pot + pot_local
1030		!! we assume the c code will do the allreduce to get the total potential
1031		!! we could do it right here if we needed to...
1032		endif
1033	<
1033	>
1034		if (do_stress) then
1035		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1036		mpi_comm_world,mpi_err)
1037		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1038		mpi_comm_world,mpi_err)
1039		endif
1040	<
1040	>
1041		#else
1042	<
1042	>
1043		if (do_stress) then
1044		tau = tau_Temp
1045		virial = virial_Temp
1046		endif
1047	<
1047	>
1048		#endif
1049	<
1049	>
1050		end subroutine do_force_loop
1051	<
1051	>
1052		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1053	<	u_l, A, f, t, pot, vpair, fpair)
1053	>	eFrame, A, f, t, pot, vpair, fpair)
1054
1055		real( kind = dp ) :: pot, vpair, sw
1056		real( kind = dp ), dimension(3) :: fpair
1057		real( kind = dp ), dimension(nLocal)   :: mfact
1058	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1058	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1059		real( kind = dp ), dimension(9,nLocal) :: A
1060		real( kind = dp ), dimension(3,nLocal) :: f
1061		real( kind = dp ), dimension(3,nLocal) :: t
#	Line 826 | Line 1065 | contains
1065		real ( kind = dp ), intent(inout) :: rijsq
1066		real ( kind = dp )                :: r
1067		real ( kind = dp ), intent(inout) :: d(3)
1068	+	real ( kind = dp ) :: ebalance
1069		integer :: me_i, me_j
1070
1071	+	integer :: iMap
1072	+
1073		r = sqrt(rijsq)
1074		vpair = 0.0d0
1075		fpair(1:3) = 0.0d0
#	Line 839 | Line 1081 | contains
1081		me_i = atid(i)
1082		me_j = atid(j)
1083		#endif
842	–
843	–	if (FF_uses_LJ .and. SIM_uses_LJ) then
844	–
845	–	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
846	–	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
847	–	endif
848	–
849	–	endif
850	–
851	–	if (FF_uses_charges .and. SIM_uses_charges) then
852	–
853	–	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
854	–	call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
855	–	endif
856	–
857	–	endif
858	–
859	–	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
860	–
861	–	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
862	–	call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
863	–	do_pot)
864	–	if (FF_uses_RF .and. SIM_uses_RF) then
865	–	call accumulate_rf(i, j, r, u_l, sw)
866	–	call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
867	–	endif
868	–	endif
1084
1085	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
1086	+
1087	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1088	+	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1089		endif
1090
1091	<	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1091	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1092	>	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093	>	pot, eFrame, f, t, do_pot)
1094
1095	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1096	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
1097	<	do_pot)
1095	>	if (FF_uses_RF .and. SIM_uses_RF) then
1096	>
1097	>	! CHECK ME (RF needs to know about all electrostatic types)
1098	>	call accumulate_rf(i, j, r, eFrame, sw)
1099	>	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100		endif
1101
1102		endif
1103
1104	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1105	+	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106	+	pot, A, f, t, do_pot)
1107	+	endif
1108
1109	<	if (FF_uses_GB .and. SIM_uses_GB) then
1110	<
1111	<	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1112	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
886	<	do_pot)
887	<	endif
1109	>	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1110	>	call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1111	>	pot, A, f, t, do_pot)
1112	>	endif
1113
1114	+	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1115	+	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1116	+	pot, A, f, t, do_pot)
1117		endif
1118	<
1119	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1120	<
1121	<	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
894	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
895	<	do_pot)
896	<	endif
897	<
1118	>
1119	>	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1120	>	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1121	>	!           pot, A, f, t, do_pot)
1122		endif
1123	+
1124	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1125	+	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1126	+	do_pot)
1127	+	endif
1128	+
1129	+	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1130	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1131	+	pot, A, f, t, do_pot)
1132	+	endif
1133	+
1134	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1135	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1136	+	pot, A, f, t, do_pot)
1137	+	endif
1138
1139		end subroutine do_pair
1140
1141		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1142	<	do_pot, do_stress, u_l, A, f, t, pot)
1142	>	do_pot, do_stress, eFrame, A, f, t, pot)
1143
1144	<	real( kind = dp ) :: pot, sw
1145	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1146	<	real (kind=dp), dimension(9,nLocal) :: A
1147	<	real (kind=dp), dimension(3,nLocal) :: f
1148	<	real (kind=dp), dimension(3,nLocal) :: t
910	<
911	<	logical, intent(inout) :: do_pot, do_stress
912	<	integer, intent(in) :: i, j
913	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
914	<	real ( kind = dp )                :: r, rc
915	<	real ( kind = dp ), intent(inout) :: d(3), dc(3)
916	<
917	<	logical :: is_EAM_i, is_EAM_j
918	<
919	<	integer :: me_i, me_j
920	<
1144	>	real( kind = dp ) :: pot, sw
1145	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1146	>	real (kind=dp), dimension(9,nLocal) :: A
1147	>	real (kind=dp), dimension(3,nLocal) :: f
1148	>	real (kind=dp), dimension(3,nLocal) :: t
1149
1150	<	r = sqrt(rijsq)
1151	<	if (SIM_uses_molecular_cutoffs) then
1152	<	rc = sqrt(rcijsq)
1153	<	else
1154	<	rc = r
927	<	endif
928	<
1150	>	logical, intent(inout) :: do_pot, do_stress
1151	>	integer, intent(in) :: i, j
1152	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1153	>	real ( kind = dp )                :: r, rc
1154	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1155
1156	+	integer :: me_i, me_j, iMap
1157	+
1158		#ifdef IS_MPI
1159	<	me_i = atid_row(i)
1160	<	me_j = atid_col(j)
1159	>	me_i = atid_row(i)
1160	>	me_j = atid_col(j)
1161		#else
1162	<	me_i = atid(i)
1163	<	me_j = atid(j)
1162	>	me_i = atid(i)
1163	>	me_j = atid(j)
1164		#endif
1165	<
1166	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1167	<
1168	<	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1169	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1170	<
1171	<	endif
1172	<
1173	<	end subroutine do_prepair
1174	<
1175	<
1176	<	subroutine do_preforce(nlocal,pot)
1177	<	integer :: nlocal
1178	<	real( kind = dp ) :: pot
1179	<
1180	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1181	<	call calc_EAM_preforce_Frho(nlocal,pot)
1182	<	endif
1183	<
1184	<
1185	<	end subroutine do_preforce
1186	<
1187	<
1188	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1189	<
1190	<	real (kind = dp), dimension(3) :: q_i
1191	<	real (kind = dp), dimension(3) :: q_j
1192	<	real ( kind = dp ), intent(out) :: r_sq
1193	<	real( kind = dp ) :: d(3), scaled(3)
1194	<	integer i
1195	<
1196	<	d(1:3) = q_j(1:3) - q_i(1:3)
1197	<
1198	<	! Wrap back into periodic box if necessary
1199	<	if ( SIM_uses_PBC ) then
1200	<
1201	<	if( .not.boxIsOrthorhombic ) then
1202	<	! calc the scaled coordinates.
1203	<
1204	<	scaled = matmul(HmatInv, d)
1205	<
1206	<	! wrap the scaled coordinates
1207	<
1208	<	scaled = scaled  - anint(scaled)
1209	<
1210	<
1211	<	! calc the wrapped real coordinates from the wrapped scaled
1212	<	! coordinates
1213	<
1214	<	d = matmul(Hmat,scaled)
1215	<
1216	<	else
1217	<	! calc the scaled coordinates.
1218	<
1219	<	do i = 1, 3
1220	<	scaled(i) = d(i) * HmatInv(i,i)
1221	<
1222	<	! wrap the scaled coordinates
1223	<
1224	<	scaled(i) = scaled(i) - anint(scaled(i))
1225	<
1226	<	! calc the wrapped real coordinates from the wrapped scaled
1227	<	! coordinates
1228	<
1229	<	d(i) = scaled(i)*Hmat(i,i)
1230	<	enddo
1231	<	endif
1232	<
1233	<	endif
1234	<
1235	<	r_sq = dot_product(d,d)
1236	<
1237	<	end subroutine get_interatomic_vector
1238	<
1239	<	subroutine zero_work_arrays()
1012	<
1165	>
1166	>	iMap = InteractionMap(me_i, me_j)%InteractionHash
1167	>
1168	>	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1169	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1170	>	endif
1171	>
1172	>	end subroutine do_prepair
1173	>
1174	>
1175	>	subroutine do_preforce(nlocal,pot)
1176	>	integer :: nlocal
1177	>	real( kind = dp ) :: pot
1178	>
1179	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1180	>	call calc_EAM_preforce_Frho(nlocal,pot)
1181	>	endif
1182	>
1183	>
1184	>	end subroutine do_preforce
1185	>
1186	>
1187	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1188	>
1189	>	real (kind = dp), dimension(3) :: q_i
1190	>	real (kind = dp), dimension(3) :: q_j
1191	>	real ( kind = dp ), intent(out) :: r_sq
1192	>	real( kind = dp ) :: d(3), scaled(3)
1193	>	integer i
1194	>
1195	>	d(1:3) = q_j(1:3) - q_i(1:3)
1196	>
1197	>	! Wrap back into periodic box if necessary
1198	>	if ( SIM_uses_PBC ) then
1199	>
1200	>	if( .not.boxIsOrthorhombic ) then
1201	>	! calc the scaled coordinates.
1202	>
1203	>	scaled = matmul(HmatInv, d)
1204	>
1205	>	! wrap the scaled coordinates
1206	>
1207	>	scaled = scaled  - anint(scaled)
1208	>
1209	>
1210	>	! calc the wrapped real coordinates from the wrapped scaled
1211	>	! coordinates
1212	>
1213	>	d = matmul(Hmat,scaled)
1214	>
1215	>	else
1216	>	! calc the scaled coordinates.
1217	>
1218	>	do i = 1, 3
1219	>	scaled(i) = d(i) * HmatInv(i,i)
1220	>
1221	>	! wrap the scaled coordinates
1222	>
1223	>	scaled(i) = scaled(i) - anint(scaled(i))
1224	>
1225	>	! calc the wrapped real coordinates from the wrapped scaled
1226	>	! coordinates
1227	>
1228	>	d(i) = scaled(i)*Hmat(i,i)
1229	>	enddo
1230	>	endif
1231	>
1232	>	endif
1233	>
1234	>	r_sq = dot_product(d,d)
1235	>
1236	>	end subroutine get_interatomic_vector
1237	>
1238	>	subroutine zero_work_arrays()
1239	>
1240		#ifdef IS_MPI
1014	–
1015	–	q_Row = 0.0_dp
1016	–	q_Col = 0.0_dp
1241
1242	<	q_group_Row = 0.0_dp
1243	<	q_group_Col = 0.0_dp
1244	<
1245	<	u_l_Row = 0.0_dp
1246	<	u_l_Col = 0.0_dp
1247	<
1248	<	A_Row = 0.0_dp
1249	<	A_Col = 0.0_dp
1250	<
1251	<	f_Row = 0.0_dp
1252	<	f_Col = 0.0_dp
1253	<	f_Temp = 0.0_dp
1254	<
1255	<	t_Row = 0.0_dp
1256	<	t_Col = 0.0_dp
1257	<	t_Temp = 0.0_dp
1258	<
1259	<	pot_Row = 0.0_dp
1260	<	pot_Col = 0.0_dp
1261	<	pot_Temp = 0.0_dp
1262	<
1263	<	rf_Row = 0.0_dp
1264	<	rf_Col = 0.0_dp
1265	<	rf_Temp = 0.0_dp
1266	<
1242	>	q_Row = 0.0_dp
1243	>	q_Col = 0.0_dp
1244	>
1245	>	q_group_Row = 0.0_dp
1246	>	q_group_Col = 0.0_dp
1247	>
1248	>	eFrame_Row = 0.0_dp
1249	>	eFrame_Col = 0.0_dp
1250	>
1251	>	A_Row = 0.0_dp
1252	>	A_Col = 0.0_dp
1253	>
1254	>	f_Row = 0.0_dp
1255	>	f_Col = 0.0_dp
1256	>	f_Temp = 0.0_dp
1257	>
1258	>	t_Row = 0.0_dp
1259	>	t_Col = 0.0_dp
1260	>	t_Temp = 0.0_dp
1261	>
1262	>	pot_Row = 0.0_dp
1263	>	pot_Col = 0.0_dp
1264	>	pot_Temp = 0.0_dp
1265	>
1266	>	rf_Row = 0.0_dp
1267	>	rf_Col = 0.0_dp
1268	>	rf_Temp = 0.0_dp
1269	>
1270		#endif
1271	<
1272	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1273	<	call clean_EAM()
1274	<	endif
1275	<
1276	<	rf = 0.0_dp
1277	<	tau_Temp = 0.0_dp
1278	<	virial_Temp = 0.0_dp
1279	<	end subroutine zero_work_arrays
1280	<
1281	<	function skipThisPair(atom1, atom2) result(skip_it)
1282	<	integer, intent(in) :: atom1
1283	<	integer, intent(in), optional :: atom2
1284	<	logical :: skip_it
1285	<	integer :: unique_id_1, unique_id_2
1286	<	integer :: me_i,me_j
1287	<	integer :: i
1288	<
1289	<	skip_it = .false.
1290	<
1291	<	!! there are a number of reasons to skip a pair or a particle
1292	<	!! mostly we do this to exclude atoms who are involved in short
1293	<	!! range interactions (bonds, bends, torsions), but we also need
1294	<	!! to exclude some overcounted interactions that result from
1295	<	!! the parallel decomposition
1296	<
1271	>
1272	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1273	>	call clean_EAM()
1274	>	endif
1275	>
1276	>	rf = 0.0_dp
1277	>	tau_Temp = 0.0_dp
1278	>	virial_Temp = 0.0_dp
1279	>	end subroutine zero_work_arrays
1280	>
1281	>	function skipThisPair(atom1, atom2) result(skip_it)
1282	>	integer, intent(in) :: atom1
1283	>	integer, intent(in), optional :: atom2
1284	>	logical :: skip_it
1285	>	integer :: unique_id_1, unique_id_2
1286	>	integer :: me_i,me_j
1287	>	integer :: i
1288	>
1289	>	skip_it = .false.
1290	>
1291	>	!! there are a number of reasons to skip a pair or a particle
1292	>	!! mostly we do this to exclude atoms who are involved in short
1293	>	!! range interactions (bonds, bends, torsions), but we also need
1294	>	!! to exclude some overcounted interactions that result from
1295	>	!! the parallel decomposition
1296	>
1297		#ifdef IS_MPI
1298	<	!! in MPI, we have to look up the unique IDs for each atom
1299	<	unique_id_1 = AtomRowToGlobal(atom1)
1298	>	!! in MPI, we have to look up the unique IDs for each atom
1299	>	unique_id_1 = AtomRowToGlobal(atom1)
1300		#else
1301	<	!! in the normal loop, the atom numbers are unique
1302	<	unique_id_1 = atom1
1301	>	!! in the normal loop, the atom numbers are unique
1302	>	unique_id_1 = atom1
1303		#endif
1304	<
1305	<	!! We were called with only one atom, so just check the global exclude
1306	<	!! list for this atom
1307	<	if (.not. present(atom2)) then
1308	<	do i = 1, nExcludes_global
1309	<	if (excludesGlobal(i) == unique_id_1) then
1310	<	skip_it = .true.
1311	<	return
1312	<	end if
1313	<	end do
1314	<	return
1315	<	end if
1316	<
1304	>
1305	>	!! We were called with only one atom, so just check the global exclude
1306	>	!! list for this atom
1307	>	if (.not. present(atom2)) then
1308	>	do i = 1, nExcludes_global
1309	>	if (excludesGlobal(i) == unique_id_1) then
1310	>	skip_it = .true.
1311	>	return
1312	>	end if
1313	>	end do
1314	>	return
1315	>	end if
1316	>
1317		#ifdef IS_MPI
1318	<	unique_id_2 = AtomColToGlobal(atom2)
1318	>	unique_id_2 = AtomColToGlobal(atom2)
1319		#else
1320	<	unique_id_2 = atom2
1320	>	unique_id_2 = atom2
1321		#endif
1322	<
1322	>
1323		#ifdef IS_MPI
1324	<	!! this situation should only arise in MPI simulations
1325	<	if (unique_id_1 == unique_id_2) then
1326	<	skip_it = .true.
1327	<	return
1328	<	end if
1329	<
1330	<	!! this prevents us from doing the pair on multiple processors
1331	<	if (unique_id_1 < unique_id_2) then
1332	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1333	<	skip_it = .true.
1334	<	return
1335	<	endif
1336	<	else
1337	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1338	<	skip_it = .true.
1339	<	return
1340	<	endif
1341	<	endif
1324	>	!! this situation should only arise in MPI simulations
1325	>	if (unique_id_1 == unique_id_2) then
1326	>	skip_it = .true.
1327	>	return
1328	>	end if
1329	>
1330	>	!! this prevents us from doing the pair on multiple processors
1331	>	if (unique_id_1 < unique_id_2) then
1332	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1333	>	skip_it = .true.
1334	>	return
1335	>	endif
1336	>	else
1337	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1338	>	skip_it = .true.
1339	>	return
1340	>	endif
1341	>	endif
1342		#endif
1343	<
1344	<	!! the rest of these situations can happen in all simulations:
1345	<	do i = 1, nExcludes_global
1346	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1347	<	(excludesGlobal(i) == unique_id_2)) then
1348	<	skip_it = .true.
1349	<	return
1350	<	endif
1351	<	enddo
1352	<
1353	<	do i = 1, nSkipsForAtom(atom1)
1354	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1355	<	skip_it = .true.
1356	<	return
1357	<	endif
1358	<	end do
1359	<
1360	<	return
1361	<	end function skipThisPair
1362	<
1363	<	function FF_UsesDirectionalAtoms() result(doesit)
1364	<	logical :: doesit
1365	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1366	<	FF_uses_GB .or. FF_uses_RF
1367	<	end function FF_UsesDirectionalAtoms
1368	<
1369	<	function FF_RequiresPrepairCalc() result(doesit)
1370	<	logical :: doesit
1371	<	doesit = FF_uses_EAM
1372	<	end function FF_RequiresPrepairCalc
1373	<
1374	<	function FF_RequiresPostpairCalc() result(doesit)
1375	<	logical :: doesit
1376	<	doesit = FF_uses_RF
1377	<	end function FF_RequiresPostpairCalc
1378	<
1343	>
1344	>	!! the rest of these situations can happen in all simulations:
1345	>	do i = 1, nExcludes_global
1346	>	if ((excludesGlobal(i) == unique_id_1) .or. &
1347	>	(excludesGlobal(i) == unique_id_2)) then
1348	>	skip_it = .true.
1349	>	return
1350	>	endif
1351	>	enddo
1352	>
1353	>	do i = 1, nSkipsForAtom(atom1)
1354	>	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1355	>	skip_it = .true.
1356	>	return
1357	>	endif
1358	>	end do
1359	>
1360	>	return
1361	>	end function skipThisPair
1362	>
1363	>	function FF_UsesDirectionalAtoms() result(doesit)
1364	>	logical :: doesit
1365	>	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366	>	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367	>	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1368	>	end function FF_UsesDirectionalAtoms
1369	>
1370	>	function FF_RequiresPrepairCalc() result(doesit)
1371	>	logical :: doesit
1372	>	doesit = FF_uses_EAM
1373	>	end function FF_RequiresPrepairCalc
1374	>
1375	>	function FF_RequiresPostpairCalc() result(doesit)
1376	>	logical :: doesit
1377	>	doesit = FF_uses_RF
1378	>	end function FF_RequiresPostpairCalc
1379	>
1380		#ifdef PROFILE
1381	<	function getforcetime() result(totalforcetime)
1382	<	real(kind=dp) :: totalforcetime
1383	<	totalforcetime = forcetime
1384	<	end function getforcetime
1381	>	function getforcetime() result(totalforcetime)
1382	>	real(kind=dp) :: totalforcetime
1383	>	totalforcetime = forcetime
1384	>	end function getforcetime
1385		#endif
1158	–
1159	–	!! This cleans componets of force arrays belonging only to fortran
1386
1387	<	subroutine add_stress_tensor(dpair, fpair)
1162	<
1163	<	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1164	<
1165	<	! because the d vector is the rj - ri vector, and
1166	<	! because fx, fy, fz are the force on atom i, we need a
1167	<	! negative sign here:
1168	<
1169	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1170	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1171	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1172	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1173	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1174	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1175	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1176	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1177	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1178	<
1179	<	virial_Temp = virial_Temp + &
1180	<	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1181	<
1182	<	end subroutine add_stress_tensor
1183	<
1184	<	end module doForces
1387	>	!! This cleans componets of force arrays belonging only to fortran
1388
1389	<	!! Interfaces for C programs to module....
1389	>	subroutine add_stress_tensor(dpair, fpair)
1390
1391	<	subroutine initFortranFF(LJMIXPOLICY, use_RF_c, thisStat)
1189	<	use doForces, ONLY: init_FF
1190	<	integer, intent(in) :: LJMIXPOLICY
1191	<	logical, intent(in) :: use_RF_c
1391	>	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1392
1393	<	integer, intent(out) :: thisStat
1394	<	call init_FF(LJMIXPOLICY, use_RF_c, thisStat)
1393	>	! because the d vector is the rj - ri vector, and
1394	>	! because fx, fy, fz are the force on atom i, we need a
1395	>	! negative sign here:
1396
1397	<	end subroutine initFortranFF
1397	>	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1398	>	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1399	>	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1400	>	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1401	>	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1402	>	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1403	>	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1404	>	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1405	>	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1406
1407	<	subroutine doForceloop(q, q_group, A, u_l, f, t, tau, pot, &
1408	<	do_pot_c, do_stress_c, error)
1200	<
1201	<	use definitions, ONLY: dp
1202	<	use simulation
1203	<	use doForces, ONLY: do_force_loop
1204	<	!! Position array provided by C, dimensioned by getNlocal
1205	<	real ( kind = dp ), dimension(3, nLocal) :: q
1206	<	!! molecular center-of-mass position array
1207	<	real ( kind = dp ), dimension(3, nGroups) :: q_group
1208	<	!! Rotation Matrix for each long range particle in simulation.
1209	<	real( kind = dp), dimension(9, nLocal) :: A
1210	<	!! Unit vectors for dipoles (lab frame)
1211	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1212	<	!! Force array provided by C, dimensioned by getNlocal
1213	<	real ( kind = dp ), dimension(3,nLocal) :: f
1214	<	!! Torsion array provided by C, dimensioned by getNlocal
1215	<	real( kind = dp ), dimension(3,nLocal) :: t
1407	>	virial_Temp = virial_Temp + &
1408	>	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1409
1410	<	!! Stress Tensor
1411	<	real( kind = dp), dimension(9) :: tau
1412	<	real ( kind = dp ) :: pot
1220	<	logical ( kind = 2) :: do_pot_c, do_stress_c
1221	<	integer :: error
1222	<
1223	<	call do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
1224	<	do_pot_c, do_stress_c, error)
1225	<
1226	<	end subroutine doForceloop
1410	>	end subroutine add_stress_tensor
1411	>
1412	>	end module doForces

Diff Legend

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