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 2266 by chuckv, Thu Jul 28 22:12:45 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.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
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	+	end subroutine createInteractionMap
251
252	<	havePropertyMap = .true.
252	>	! 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.
253	>	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
262
263	<	end subroutine createPropertyMap
263	>	stat = 0
264	>	if (.not. haveInteractionMap) then
265
266	+	call createInteractionMap(myStatus)
267	+
268	+	if (myStatus .ne. 0) then
269	+	write(default_error, *) 'createInteractionMap failed in doForces!'
270	+	stat = -1
271	+	return
272	+	endif
273	+	endif
274	+
275	+
276	+	nAtypes = getSize(atypes)
277	+	! If we pass a default rcut, set all atypes to that cutoff distance
278	+	if(present(defaultRList)) then
279	+	InteractionMap(:,:)%rList = defaultRList
280	+	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
281	+	haveRlist = .true.
282	+	return
283	+	end if
284	+
285	+	do map_i = 1,nAtypes
286	+	do map_j = map_i,nAtypes
287	+	iMap = InteractionMap(map_i, map_j)%InteractionHash
288	+
289	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
290	+	!            thisRCut = getLJCutOff(map_i,map_j)
291	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
292	+	endif
293	+
294	+	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
295	+	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
296	+	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
339	+
340	+
341	+	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
342	+	!!$  subroutine setRlistDF( this_rlist )
343	+	!!$
344	+	!!$   real(kind=dp) :: this_rlist
345	+	!!$
346	+	!!$    rlist = this_rlist
347	+	!!$    rlistsq = rlist * rlist
348	+	!!$
349	+	!!$    haveRlist = .true.
350	+	!!$
351	+	!!$  end subroutine setRlistDF
352	+
353	+
354		subroutine setSimVariables()
355	<	SIM_uses_LJ = SimUsesLJ()
356	<	SIM_uses_sticky = SimUsesSticky()
357	<	SIM_uses_charges = SimUsesCharges()
358	<	SIM_uses_dipoles = SimUsesDipoles()
359	<	SIM_uses_RF = SimUsesRF()
360	<	SIM_uses_GB = SimUsesGB()
355	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
356	>	SIM_uses_LennardJones = SimUsesLennardJones()
357	>	SIM_uses_Electrostatics = SimUsesElectrostatics()
358	>	SIM_uses_Charges = SimUsesCharges()
359	>	SIM_uses_Dipoles = SimUsesDipoles()
360	>	SIM_uses_Sticky = SimUsesSticky()
361	>	SIM_uses_StickyPower = SimUsesStickyPower()
362	>	SIM_uses_GayBerne = SimUsesGayBerne()
363		SIM_uses_EAM = SimUsesEAM()
364	+	SIM_uses_Shapes = SimUsesShapes()
365	+	SIM_uses_FLARB = SimUsesFLARB()
366	+	SIM_uses_RF = SimUsesRF()
367		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
368		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
166	–	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
369		SIM_uses_PBC = SimUsesPBC()
168	–	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
370
371		haveSIMvariables = .true.
372
#	Line 178 | Line 379 | contains
379		integer :: myStatus
380
381		error = 0
181	–
182	–	if (.not. havePropertyMap) then
382
383	+	if (.not. haveInteractionMap) then
384	+
385		myStatus = 0
386
387	<	call createPropertyMap(myStatus)
387	>	call createInteractionMap(myStatus)
388
389		if (myStatus .ne. 0) then
390	<	write(default_error, *) 'createPropertyMap failed in doForces!'
390	>	write(default_error, *) 'createInteractionMap failed in doForces!'
391		error = -1
392		return
393		endif
#	Line 202 | Line 403 | contains
403		return
404		endif
405
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	–
406		if (.not. haveNeighborList) then
407		write(default_error, *) 'neighbor list has not been initialized in doForces!'
408		error = -1
#	Line 231 | Line 424 | contains
424		#endif
425		return
426		end subroutine doReadyCheck
234	–
427
236	–	subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
428
429	<	integer, intent(in) :: LJMIXPOLICY
429	>	subroutine init_FF(use_RF_c, thisStat)
430	>
431		logical, intent(in) :: use_RF_c
432
433		integer, intent(out) :: thisStat
#	Line 248 | Line 440 | contains
440
441		!! Fortran's version of a cast:
442		FF_uses_RF = use_RF_c
443	<
443	>
444		!! init_FF is called *after* all of the atom types have been
445		!! defined in atype_module using the new_atype subroutine.
446		!!
447		!! this will scan through the known atypes and figure out what
448		!! interactions are used by the force field.
449	<
450	<	FF_uses_LJ = .false.
451	<	FF_uses_sticky = .false.
452	<	FF_uses_charges = .false.
453	<	FF_uses_dipoles = .false.
454	<	FF_uses_GB = .false.
449	>
450	>	FF_uses_DirectionalAtoms = .false.
451	>	FF_uses_LennardJones = .false.
452	>	FF_uses_Electrostatics = .false.
453	>	FF_uses_Charges = .false.
454	>	FF_uses_Dipoles = .false.
455	>	FF_uses_Sticky = .false.
456	>	FF_uses_StickyPower = .false.
457	>	FF_uses_GayBerne = .false.
458		FF_uses_EAM = .false.
459	<
460	<	call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
266	<	if (nMatches .gt. 0) FF_uses_LJ = .true.
459	>	FF_uses_Shapes = .false.
460	>	FF_uses_FLARB = .false.
461
462	<	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
463	<	if (nMatches .gt. 0) FF_uses_charges = .true.
462	>	call getMatchingElementList(atypes, "is_Directional", .true., &
463	>	nMatches, MatchList)
464	>	if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
465
466	<	call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
467	<	if (nMatches .gt. 0) FF_uses_dipoles = .true.
468	<
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	>
483	>	call getMatchingElementList(atypes, "is_Dipole", .true., &
484	>	nMatches, MatchList)
485	>	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) FF_uses_Sticky = .true.
502	<
503	<	call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
504	<	if (nMatches .gt. 0) FF_uses_GB = .true.
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
512
513	+	call getMatchingElementList(atypes, "is_GayBerne", .true., &
514	+	nMatches, MatchList)
515	+	if (nMatches .gt. 0) then
516	+	FF_uses_GayBerne = .true.
517	+	FF_uses_DirectionalAtoms = .true.
518	+	endif
519	+
520		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
521		if (nMatches .gt. 0) FF_uses_EAM = .true.
522	<
522	>
523	>	call getMatchingElementList(atypes, "is_Shape", .true., &
524	>	nMatches, MatchList)
525	>	if (nMatches .gt. 0) then
526	>	FF_uses_Shapes = .true.
527	>	FF_uses_DirectionalAtoms = .true.
528	>	endif
529	>
530	>	call getMatchingElementList(atypes, "is_FLARB", .true., &
531	>	nMatches, MatchList)
532	>	if (nMatches .gt. 0) FF_uses_FLARB = .true.
533	>
534		!! Assume sanity (for the sake of argument)
535		haveSaneForceField = .true.
536
537		!! check to make sure the FF_uses_RF setting makes sense
538	<
538	>
539		if (FF_uses_dipoles) then
540		if (FF_uses_RF) then
541		dielect = getDielect()
#	Line 298 | Line 548 | contains
548		haveSaneForceField = .false.
549		return
550		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.
551		endif
552
553	<	if (FF_uses_sticky) then
554	<	call check_sticky_FF(my_status)
555	<	if (my_status /= 0) then
556	<	thisStat = -1
557	<	haveSaneForceField = .false.
558	<	return
559	<	end if
560	<	endif
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
333	–
563		if (FF_uses_EAM) then
564	<	call init_EAM_FF(my_status)
564	>	call init_EAM_FF(my_status)
565		if (my_status /= 0) then
566		write(default_error, *) "init_EAM_FF returned a bad status"
567		thisStat = -1
#	Line 341 | Line 570 | contains
570		end if
571		endif
572
573	<	if (FF_uses_GB) then
573	>	if (FF_uses_GayBerne) then
574		call check_gb_pair_FF(my_status)
575		if (my_status .ne. 0) then
576		thisStat = -1
#	Line 350 | Line 579 | contains
579		endif
580		endif
581
582	<	if (FF_uses_GB .and. FF_uses_LJ) then
582	>	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583		endif
584	+
585		if (.not. haveNeighborList) then
586		!! Create neighbor lists
587		call expandNeighborList(nLocal, my_status)
#	Line 363 | Line 593 | contains
593		haveNeighborList = .true.
594		endif
595
366	–
367	–
596		end subroutine init_FF
369	–
597
598	+
599		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
600		!------------------------------------------------------------->
601	<	subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
601	>	subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
602		do_pot_c, do_stress_c, error)
603		!! Position array provided by C, dimensioned by getNlocal
604		real ( kind = dp ), dimension(3, nLocal) :: q
#	Line 379 | Line 607 | contains
607		!! Rotation Matrix for each long range particle in simulation.
608		real( kind = dp), dimension(9, nLocal) :: A
609		!! Unit vectors for dipoles (lab frame)
610	<	real( kind = dp ), dimension(3,nLocal) :: u_l
610	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
611		!! Force array provided by C, dimensioned by getNlocal
612		real ( kind = dp ), dimension(3,nLocal) :: f
613		!! Torsion array provided by C, dimensioned by getNlocal
#	Line 417 | Line 645 | contains
645		integer :: localError
646		integer :: propPack_i, propPack_j
647		integer :: loopStart, loopEnd, loop
648	<
648	>	integer :: iMap
649		real(kind=dp) :: listSkin = 1.0
650	<
650	>
651		!! initialize local variables
652	<
652	>
653		#ifdef IS_MPI
654		pot_local = 0.0_dp
655		nAtomsInRow   = getNatomsInRow(plan_atom_row)
#	Line 431 | Line 659 | contains
659		#else
660		natoms = nlocal
661		#endif
662	<
662	>
663		call doReadyCheck(localError)
664		if ( localError .ne. 0 ) then
665		call handleError("do_force_loop", "Not Initialized")
#	Line 439 | Line 667 | contains
667		return
668		end if
669		call zero_work_arrays()
670	<
670	>
671		do_pot = do_pot_c
672		do_stress = do_stress_c
673	<
673	>
674		! Gather all information needed by all force loops:
675	<
675	>
676		#ifdef IS_MPI
677	<
677	>
678		call gather(q, q_Row, plan_atom_row_3d)
679		call gather(q, q_Col, plan_atom_col_3d)
680
681		call gather(q_group, q_group_Row, plan_group_row_3d)
682		call gather(q_group, q_group_Col, plan_group_col_3d)
683	<
684	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
685	<	call gather(u_l, u_l_Row, plan_atom_row_3d)
686	<	call gather(u_l, u_l_Col, plan_atom_col_3d)
687	<
683	>
684	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
685	>	call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
686	>	call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
687	>
688		call gather(A, A_Row, plan_atom_row_rotation)
689		call gather(A, A_Col, plan_atom_col_rotation)
690		endif
691	<
691	>
692		#endif
693	<
693	>
694		!! Begin force loop timing:
695		#ifdef PROFILE
696		call cpu_time(forceTimeInitial)
697		nloops = nloops + 1
698		#endif
699	<
699	>
700		loopEnd = PAIR_LOOP
701		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
702		loopStart = PREPAIR_LOOP
#	Line 483 | Line 711 | contains
711		if (loop .eq. loopStart) then
712		#ifdef IS_MPI
713		call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
714	<	update_nlist)
714	>	update_nlist)
715		#else
716		call checkNeighborList(nGroups, q_group, listSkin, &
717	<	update_nlist)
717	>	update_nlist)
718		#endif
719		endif
720	<
720	>
721		if (update_nlist) then
722		!! save current configuration and construct neighbor list
723		#ifdef IS_MPI
#	Line 500 | Line 728 | contains
728		neighborListSize = size(list)
729		nlist = 0
730		endif
731	<
731	>
732		istart = 1
733		#ifdef IS_MPI
734		iend = nGroupsInRow
#	Line 510 | Line 738 | contains
738		outer: do i = istart, iend
739
740		if (update_nlist) point(i) = nlist + 1
741	<
741	>
742		n_in_i = groupStartRow(i+1) - groupStartRow(i)
743	<
743	>
744		if (update_nlist) then
745		#ifdef IS_MPI
746	+	me_i = atid_row(i)
747		jstart = 1
748		jend = nGroupsInCol
749		#else
750	+	me_i = atid(i)
751		jstart = i+1
752		jend = nGroups
753		#endif
#	Line 527 | Line 757 | contains
757		! make sure group i has neighbors
758		if (jstart .gt. jend) cycle outer
759		endif
760	<
760	>
761		do jnab = jstart, jend
762		if (update_nlist) then
763		j = jnab
#	Line 536 | Line 766 | contains
766		endif
767
768		#ifdef IS_MPI
769	+	me_j = atid_col(j)
770		call get_interatomic_vector(q_group_Row(:,i), &
771		q_group_Col(:,j), d_grp, rgrpsq)
772		#else
773	+	me_j = atid(j)
774		call get_interatomic_vector(q_group(:,i), &
775		q_group(:,j), d_grp, rgrpsq)
776		#endif
777
778	<	if (rgrpsq < rlistsq) then
778	>	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
779		if (update_nlist) then
780		nlist = nlist + 1
781	<
781	>
782		if (nlist > neighborListSize) then
783		#ifdef IS_MPI
784		call expandNeighborList(nGroupsInRow, listerror)
#	Line 560 | Line 792 | contains
792		end if
793		neighborListSize = size(list)
794		endif
795	<
795	>
796		list(nlist) = j
797		endif
798	<
798	>
799		if (loop .eq. PAIR_LOOP) then
800		vij = 0.0d0
801		fij(1:3) = 0.0d0
802		endif
803	<
803	>
804		call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
805		in_switching_region)
806	<
806	>
807		n_in_j = groupStartCol(j+1) - groupStartCol(j)
808	<
808	>
809		do ia = groupStartRow(i), groupStartRow(i+1)-1
810	<
810	>
811		atom1 = groupListRow(ia)
812	<
812	>
813		inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
814	<
814	>
815		atom2 = groupListCol(jb)
816	<
816	>
817		if (skipThisPair(atom1, atom2)) cycle inner
818
819		if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
#	Line 601 | Line 833 | contains
833		#ifdef IS_MPI
834		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
835		rgrpsq, d_grp, do_pot, do_stress, &
836	<	u_l, A, f, t, pot_local)
836	>	eFrame, A, f, t, pot_local)
837		#else
838		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
839		rgrpsq, d_grp, do_pot, do_stress, &
840	<	u_l, A, f, t, pot)
840	>	eFrame, A, f, t, pot)
841		#endif
842		else
843		#ifdef IS_MPI
844		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
845		do_pot, &
846	<	u_l, A, f, t, pot_local, vpair, fpair)
846	>	eFrame, A, f, t, pot_local, vpair, fpair)
847		#else
848		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
849		do_pot,  &
850	<	u_l, A, f, t, pot, vpair, fpair)
850	>	eFrame, A, f, t, pot, vpair, fpair)
851		#endif
852
853		vij = vij + vpair
#	Line 623 | Line 855 | contains
855		endif
856		enddo inner
857		enddo
858	<
858	>
859		if (loop .eq. PAIR_LOOP) then
860		if (in_switching_region) then
861		swderiv = vij*dswdr/rgrp
862		fij(1) = fij(1) + swderiv*d_grp(1)
863		fij(2) = fij(2) + swderiv*d_grp(2)
864		fij(3) = fij(3) + swderiv*d_grp(3)
865	<
865	>
866		do ia=groupStartRow(i), groupStartRow(i+1)-1
867		atom1=groupListRow(ia)
868		mf = mfactRow(atom1)
#	Line 644 | Line 876 | contains
876		f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
877		#endif
878		enddo
879	<
879	>
880		do jb=groupStartCol(j), groupStartCol(j+1)-1
881		atom2=groupListCol(jb)
882		mf = mfactCol(atom2)
#	Line 659 | Line 891 | contains
891		#endif
892		enddo
893		endif
894	<
894	>
895		if (do_stress) call add_stress_tensor(d_grp, fij)
896		endif
897		end if
898		enddo
899		enddo outer
900	<
900	>
901		if (update_nlist) then
902		#ifdef IS_MPI
903		point(nGroupsInRow + 1) = nlist + 1
#	Line 679 | Line 911 | contains
911		update_nlist = .false.
912		endif
913		endif
914	<
914	>
915		if (loop .eq. PREPAIR_LOOP) then
916		call do_preforce(nlocal, pot)
917		endif
918	<
918	>
919		enddo
920	<
920	>
921		!! Do timing
922		#ifdef PROFILE
923		call cpu_time(forceTimeFinal)
924		forceTime = forceTime + forceTimeFinal - forceTimeInitial
925		#endif
926	<
926	>
927		#ifdef IS_MPI
928		!!distribute forces
929	<
929	>
930		f_temp = 0.0_dp
931		call scatter(f_Row,f_temp,plan_atom_row_3d)
932		do i = 1,nlocal
933		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
934		end do
935	<
935	>
936		f_temp = 0.0_dp
937		call scatter(f_Col,f_temp,plan_atom_col_3d)
938		do i = 1,nlocal
939		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
940		end do
941	<
942	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
941	>
942	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
943		t_temp = 0.0_dp
944		call scatter(t_Row,t_temp,plan_atom_row_3d)
945		do i = 1,nlocal
#	Line 715 | Line 947 | contains
947		end do
948		t_temp = 0.0_dp
949		call scatter(t_Col,t_temp,plan_atom_col_3d)
950	<
950	>
951		do i = 1,nlocal
952		t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
953		end do
954		endif
955	<
955	>
956		if (do_pot) then
957		! scatter/gather pot_row into the members of my column
958		call scatter(pot_Row, pot_Temp, plan_atom_row)
959	<
959	>
960		! scatter/gather pot_local into all other procs
961		! add resultant to get total pot
962		do i = 1, nlocal
963		pot_local = pot_local + pot_Temp(i)
964		enddo
965	<
965	>
966		pot_Temp = 0.0_DP
967	<
967	>
968		call scatter(pot_Col, pot_Temp, plan_atom_col)
969		do i = 1, nlocal
970		pot_local = pot_local + pot_Temp(i)
971		enddo
972	<
972	>
973		endif
974		#endif
975	<
975	>
976		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
977	<
977	>
978		if (FF_uses_RF .and. SIM_uses_RF) then
979	<
979	>
980		#ifdef IS_MPI
981		call scatter(rf_Row,rf,plan_atom_row_3d)
982		call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
#	Line 752 | Line 984 | contains
984		rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
985		end do
986		#endif
987	<
987	>
988		do i = 1, nLocal
989	<
989	>
990		rfpot = 0.0_DP
991		#ifdef IS_MPI
992		me_i = atid_row(i)
993		#else
994		me_i = atid(i)
995		#endif
996	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
997
998	<	if (PropertyMap(me_i)%is_DP) then
999	<
1000	<	mu_i = PropertyMap(me_i)%dipole_moment
1001	<
998	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
999	>
1000	>	mu_i = getDipoleMoment(me_i)
1001	>
1002		!! The reaction field needs to include a self contribution
1003		!! to the field:
1004	<	call accumulate_self_rf(i, mu_i, u_l)
1004	>	call accumulate_self_rf(i, mu_i, eFrame)
1005		!! Get the reaction field contribution to the
1006		!! potential and torques:
1007	<	call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
1007	>	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1008		#ifdef IS_MPI
1009		pot_local = pot_local + rfpot
1010		#else
1011		pot = pot + rfpot
1012	<
1012	>
1013		#endif
1014	<	endif
1014	>	endif
1015		enddo
1016		endif
1017		endif
1018	<
1019	<
1018	>
1019	>
1020		#ifdef IS_MPI
1021	<
1021	>
1022		if (do_pot) then
1023		pot = pot + pot_local
1024		!! we assume the c code will do the allreduce to get the total potential
1025		!! we could do it right here if we needed to...
1026		endif
1027	<
1027	>
1028		if (do_stress) then
1029		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1030		mpi_comm_world,mpi_err)
1031		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1032		mpi_comm_world,mpi_err)
1033		endif
1034	<
1034	>
1035		#else
1036	<
1036	>
1037		if (do_stress) then
1038		tau = tau_Temp
1039		virial = virial_Temp
1040		endif
1041	<
1041	>
1042		#endif
1043	<
1043	>
1044		end subroutine do_force_loop
1045	<
1045	>
1046		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1047	<	u_l, A, f, t, pot, vpair, fpair)
1047	>	eFrame, A, f, t, pot, vpair, fpair)
1048
1049		real( kind = dp ) :: pot, vpair, sw
1050		real( kind = dp ), dimension(3) :: fpair
1051		real( kind = dp ), dimension(nLocal)   :: mfact
1052	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1052	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1053		real( kind = dp ), dimension(9,nLocal) :: A
1054		real( kind = dp ), dimension(3,nLocal) :: f
1055		real( kind = dp ), dimension(3,nLocal) :: t
#	Line 826 | Line 1059 | contains
1059		real ( kind = dp ), intent(inout) :: rijsq
1060		real ( kind = dp )                :: r
1061		real ( kind = dp ), intent(inout) :: d(3)
1062	+	real ( kind = dp ) :: ebalance
1063		integer :: me_i, me_j
1064
1065	+	integer :: iMap
1066	+
1067		r = sqrt(rijsq)
1068		vpair = 0.0d0
1069		fpair(1:3) = 0.0d0
#	Line 839 | Line 1075 | contains
1075		me_i = atid(i)
1076		me_j = atid(j)
1077		#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
1078
1079	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
1080	+
1081	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1082	+	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1083		endif
1084
1085	<	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1085	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1086	>	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087	>	pot, eFrame, f, t, do_pot)
1088
1089	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1090	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
1091	<	do_pot)
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
1096		endif
1097
1098	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1099	+	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1100	+	pot, A, f, t, do_pot)
1101	+	endif
1102
1103	<	if (FF_uses_GB .and. SIM_uses_GB) then
1104	<
1105	<	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1106	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
886	<	do_pot)
887	<	endif
1103	>	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1104	>	call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1105	>	pot, A, f, t, do_pot)
1106	>	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
1112	<
1113	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1114	<
1115	<	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	<
1112	>
1113	>	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1114	>	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1115	>	!           pot, A, f, t, do_pot)
1116		endif
1117	+
1118	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1119	+	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1120	+	do_pot)
1121	+	endif
1122	+
1123	+	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1124	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1125	+	pot, A, f, t, do_pot)
1126	+	endif
1127	+
1128	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1129	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1130	+	pot, A, f, t, do_pot)
1131	+	endif
1132
1133		end subroutine do_pair
1134
1135		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1136	<	do_pot, do_stress, u_l, A, f, t, pot)
1136	>	do_pot, do_stress, eFrame, A, f, t, pot)
1137
1138	<	real( kind = dp ) :: pot, sw
1139	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1140	<	real (kind=dp), dimension(9,nLocal) :: A
1141	<	real (kind=dp), dimension(3,nLocal) :: f
1142	<	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	<
1138	>	real( kind = dp ) :: pot, sw
1139	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1140	>	real (kind=dp), dimension(9,nLocal) :: A
1141	>	real (kind=dp), dimension(3,nLocal) :: f
1142	>	real (kind=dp), dimension(3,nLocal) :: t
1143
1144	<	r = sqrt(rijsq)
1145	<	if (SIM_uses_molecular_cutoffs) then
1146	<	rc = sqrt(rcijsq)
1147	<	else
1148	<	rc = r
927	<	endif
928	<
1144	>	logical, intent(inout) :: do_pot, do_stress
1145	>	integer, intent(in) :: i, j
1146	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1147	>	real ( kind = dp )                :: r, rc
1148	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1149
1150	+	integer :: me_i, me_j, iMap
1151	+
1152		#ifdef IS_MPI
1153	<	me_i = atid_row(i)
1154	<	me_j = atid_col(j)
1153	>	me_i = atid_row(i)
1154	>	me_j = atid_col(j)
1155		#else
1156	<	me_i = atid(i)
1157	<	me_j = atid(j)
1156	>	me_i = atid(i)
1157	>	me_j = atid(j)
1158		#endif
1159	<
1160	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1161	<
1162	<	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1163	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1164	<
1165	<	endif
1166	<
1167	<	end subroutine do_prepair
1168	<
1169	<
1170	<	subroutine do_preforce(nlocal,pot)
1171	<	integer :: nlocal
1172	<	real( kind = dp ) :: pot
1173	<
1174	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1175	<	call calc_EAM_preforce_Frho(nlocal,pot)
1176	<	endif
1177	<
1178	<
1179	<	end subroutine do_preforce
1180	<
1181	<
1182	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1183	<
1184	<	real (kind = dp), dimension(3) :: q_i
1185	<	real (kind = dp), dimension(3) :: q_j
1186	<	real ( kind = dp ), intent(out) :: r_sq
1187	<	real( kind = dp ) :: d(3), scaled(3)
1188	<	integer i
1189	<
1190	<	d(1:3) = q_j(1:3) - q_i(1:3)
1191	<
1192	<	! Wrap back into periodic box if necessary
1193	<	if ( SIM_uses_PBC ) then
1194	<
1195	<	if( .not.boxIsOrthorhombic ) then
1196	<	! calc the scaled coordinates.
1197	<
1198	<	scaled = matmul(HmatInv, d)
1199	<
1200	<	! wrap the scaled coordinates
1201	<
1202	<	scaled = scaled  - anint(scaled)
1203	<
1204	<
1205	<	! calc the wrapped real coordinates from the wrapped scaled
1206	<	! coordinates
1207	<
1208	<	d = matmul(Hmat,scaled)
1209	<
1210	<	else
1211	<	! calc the scaled coordinates.
1212	<
1213	<	do i = 1, 3
1214	<	scaled(i) = d(i) * HmatInv(i,i)
1215	<
1216	<	! wrap the scaled coordinates
1217	<
1218	<	scaled(i) = scaled(i) - anint(scaled(i))
1219	<
1220	<	! calc the wrapped real coordinates from the wrapped scaled
1221	<	! coordinates
1222	<
1223	<	d(i) = scaled(i)*Hmat(i,i)
1224	<	enddo
1225	<	endif
1226	<
1227	<	endif
1228	<
1229	<	r_sq = dot_product(d,d)
1230	<
1231	<	end subroutine get_interatomic_vector
1232	<
1233	<	subroutine zero_work_arrays()
1012	<
1159	>
1160	>	iMap = InteractionMap(me_i, me_j)%InteractionHash
1161	>
1162	>	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1163	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1164	>	endif
1165	>
1166	>	end subroutine do_prepair
1167	>
1168	>
1169	>	subroutine do_preforce(nlocal,pot)
1170	>	integer :: nlocal
1171	>	real( kind = dp ) :: pot
1172	>
1173	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1174	>	call calc_EAM_preforce_Frho(nlocal,pot)
1175	>	endif
1176	>
1177	>
1178	>	end subroutine do_preforce
1179	>
1180	>
1181	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1182	>
1183	>	real (kind = dp), dimension(3) :: q_i
1184	>	real (kind = dp), dimension(3) :: q_j
1185	>	real ( kind = dp ), intent(out) :: r_sq
1186	>	real( kind = dp ) :: d(3), scaled(3)
1187	>	integer i
1188	>
1189	>	d(1:3) = q_j(1:3) - q_i(1:3)
1190	>
1191	>	! Wrap back into periodic box if necessary
1192	>	if ( SIM_uses_PBC ) then
1193	>
1194	>	if( .not.boxIsOrthorhombic ) then
1195	>	! calc the scaled coordinates.
1196	>
1197	>	scaled = matmul(HmatInv, d)
1198	>
1199	>	! wrap the scaled coordinates
1200	>
1201	>	scaled = scaled  - anint(scaled)
1202	>
1203	>
1204	>	! calc the wrapped real coordinates from the wrapped scaled
1205	>	! coordinates
1206	>
1207	>	d = matmul(Hmat,scaled)
1208	>
1209	>	else
1210	>	! calc the scaled coordinates.
1211	>
1212	>	do i = 1, 3
1213	>	scaled(i) = d(i) * HmatInv(i,i)
1214	>
1215	>	! wrap the scaled coordinates
1216	>
1217	>	scaled(i) = scaled(i) - anint(scaled(i))
1218	>
1219	>	! calc the wrapped real coordinates from the wrapped scaled
1220	>	! coordinates
1221	>
1222	>	d(i) = scaled(i)*Hmat(i,i)
1223	>	enddo
1224	>	endif
1225	>
1226	>	endif
1227	>
1228	>	r_sq = dot_product(d,d)
1229	>
1230	>	end subroutine get_interatomic_vector
1231	>
1232	>	subroutine zero_work_arrays()
1233	>
1234		#ifdef IS_MPI
1014	–
1015	–	q_Row = 0.0_dp
1016	–	q_Col = 0.0_dp
1235
1236	<	q_group_Row = 0.0_dp
1237	<	q_group_Col = 0.0_dp
1238	<
1239	<	u_l_Row = 0.0_dp
1240	<	u_l_Col = 0.0_dp
1241	<
1242	<	A_Row = 0.0_dp
1243	<	A_Col = 0.0_dp
1244	<
1245	<	f_Row = 0.0_dp
1246	<	f_Col = 0.0_dp
1247	<	f_Temp = 0.0_dp
1248	<
1249	<	t_Row = 0.0_dp
1250	<	t_Col = 0.0_dp
1251	<	t_Temp = 0.0_dp
1252	<
1253	<	pot_Row = 0.0_dp
1254	<	pot_Col = 0.0_dp
1255	<	pot_Temp = 0.0_dp
1256	<
1257	<	rf_Row = 0.0_dp
1258	<	rf_Col = 0.0_dp
1259	<	rf_Temp = 0.0_dp
1260	<
1236	>	q_Row = 0.0_dp
1237	>	q_Col = 0.0_dp
1238	>
1239	>	q_group_Row = 0.0_dp
1240	>	q_group_Col = 0.0_dp
1241	>
1242	>	eFrame_Row = 0.0_dp
1243	>	eFrame_Col = 0.0_dp
1244	>
1245	>	A_Row = 0.0_dp
1246	>	A_Col = 0.0_dp
1247	>
1248	>	f_Row = 0.0_dp
1249	>	f_Col = 0.0_dp
1250	>	f_Temp = 0.0_dp
1251	>
1252	>	t_Row = 0.0_dp
1253	>	t_Col = 0.0_dp
1254	>	t_Temp = 0.0_dp
1255	>
1256	>	pot_Row = 0.0_dp
1257	>	pot_Col = 0.0_dp
1258	>	pot_Temp = 0.0_dp
1259	>
1260	>	rf_Row = 0.0_dp
1261	>	rf_Col = 0.0_dp
1262	>	rf_Temp = 0.0_dp
1263	>
1264		#endif
1265	<
1266	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1267	<	call clean_EAM()
1268	<	endif
1269	<
1270	<	rf = 0.0_dp
1271	<	tau_Temp = 0.0_dp
1272	<	virial_Temp = 0.0_dp
1273	<	end subroutine zero_work_arrays
1274	<
1275	<	function skipThisPair(atom1, atom2) result(skip_it)
1276	<	integer, intent(in) :: atom1
1277	<	integer, intent(in), optional :: atom2
1278	<	logical :: skip_it
1279	<	integer :: unique_id_1, unique_id_2
1280	<	integer :: me_i,me_j
1281	<	integer :: i
1282	<
1283	<	skip_it = .false.
1284	<
1285	<	!! there are a number of reasons to skip a pair or a particle
1286	<	!! mostly we do this to exclude atoms who are involved in short
1287	<	!! range interactions (bonds, bends, torsions), but we also need
1288	<	!! to exclude some overcounted interactions that result from
1289	<	!! the parallel decomposition
1290	<
1265	>
1266	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1267	>	call clean_EAM()
1268	>	endif
1269	>
1270	>	rf = 0.0_dp
1271	>	tau_Temp = 0.0_dp
1272	>	virial_Temp = 0.0_dp
1273	>	end subroutine zero_work_arrays
1274	>
1275	>	function skipThisPair(atom1, atom2) result(skip_it)
1276	>	integer, intent(in) :: atom1
1277	>	integer, intent(in), optional :: atom2
1278	>	logical :: skip_it
1279	>	integer :: unique_id_1, unique_id_2
1280	>	integer :: me_i,me_j
1281	>	integer :: i
1282	>
1283	>	skip_it = .false.
1284	>
1285	>	!! there are a number of reasons to skip a pair or a particle
1286	>	!! mostly we do this to exclude atoms who are involved in short
1287	>	!! range interactions (bonds, bends, torsions), but we also need
1288	>	!! to exclude some overcounted interactions that result from
1289	>	!! the parallel decomposition
1290	>
1291		#ifdef IS_MPI
1292	<	!! in MPI, we have to look up the unique IDs for each atom
1293	<	unique_id_1 = AtomRowToGlobal(atom1)
1292	>	!! in MPI, we have to look up the unique IDs for each atom
1293	>	unique_id_1 = AtomRowToGlobal(atom1)
1294		#else
1295	<	!! in the normal loop, the atom numbers are unique
1296	<	unique_id_1 = atom1
1295	>	!! in the normal loop, the atom numbers are unique
1296	>	unique_id_1 = atom1
1297		#endif
1298	<
1299	<	!! We were called with only one atom, so just check the global exclude
1300	<	!! list for this atom
1301	<	if (.not. present(atom2)) then
1302	<	do i = 1, nExcludes_global
1303	<	if (excludesGlobal(i) == unique_id_1) then
1304	<	skip_it = .true.
1305	<	return
1306	<	end if
1307	<	end do
1308	<	return
1309	<	end if
1310	<
1298	>
1299	>	!! We were called with only one atom, so just check the global exclude
1300	>	!! list for this atom
1301	>	if (.not. present(atom2)) then
1302	>	do i = 1, nExcludes_global
1303	>	if (excludesGlobal(i) == unique_id_1) then
1304	>	skip_it = .true.
1305	>	return
1306	>	end if
1307	>	end do
1308	>	return
1309	>	end if
1310	>
1311		#ifdef IS_MPI
1312	<	unique_id_2 = AtomColToGlobal(atom2)
1312	>	unique_id_2 = AtomColToGlobal(atom2)
1313		#else
1314	<	unique_id_2 = atom2
1314	>	unique_id_2 = atom2
1315		#endif
1316	<
1316	>
1317		#ifdef IS_MPI
1318	<	!! this situation should only arise in MPI simulations
1319	<	if (unique_id_1 == unique_id_2) then
1320	<	skip_it = .true.
1321	<	return
1322	<	end if
1323	<
1324	<	!! this prevents us from doing the pair on multiple processors
1325	<	if (unique_id_1 < unique_id_2) then
1326	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1327	<	skip_it = .true.
1328	<	return
1329	<	endif
1330	<	else
1331	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1332	<	skip_it = .true.
1333	<	return
1334	<	endif
1335	<	endif
1318	>	!! this situation should only arise in MPI simulations
1319	>	if (unique_id_1 == unique_id_2) then
1320	>	skip_it = .true.
1321	>	return
1322	>	end if
1323	>
1324	>	!! this prevents us from doing the pair on multiple processors
1325	>	if (unique_id_1 < unique_id_2) then
1326	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1327	>	skip_it = .true.
1328	>	return
1329	>	endif
1330	>	else
1331	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1332	>	skip_it = .true.
1333	>	return
1334	>	endif
1335	>	endif
1336		#endif
1337	<
1338	<	!! the rest of these situations can happen in all simulations:
1339	<	do i = 1, nExcludes_global
1340	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1341	<	(excludesGlobal(i) == unique_id_2)) then
1342	<	skip_it = .true.
1343	<	return
1344	<	endif
1345	<	enddo
1346	<
1347	<	do i = 1, nSkipsForAtom(atom1)
1348	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1349	<	skip_it = .true.
1350	<	return
1351	<	endif
1352	<	end do
1353	<
1354	<	return
1355	<	end function skipThisPair
1356	<
1357	<	function FF_UsesDirectionalAtoms() result(doesit)
1358	<	logical :: doesit
1359	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1360	<	FF_uses_GB .or. FF_uses_RF
1361	<	end function FF_UsesDirectionalAtoms
1362	<
1363	<	function FF_RequiresPrepairCalc() result(doesit)
1364	<	logical :: doesit
1365	<	doesit = FF_uses_EAM
1366	<	end function FF_RequiresPrepairCalc
1367	<
1368	<	function FF_RequiresPostpairCalc() result(doesit)
1369	<	logical :: doesit
1370	<	doesit = FF_uses_RF
1371	<	end function FF_RequiresPostpairCalc
1372	<
1337	>
1338	>	!! the rest of these situations can happen in all simulations:
1339	>	do i = 1, nExcludes_global
1340	>	if ((excludesGlobal(i) == unique_id_1) .or. &
1341	>	(excludesGlobal(i) == unique_id_2)) then
1342	>	skip_it = .true.
1343	>	return
1344	>	endif
1345	>	enddo
1346	>
1347	>	do i = 1, nSkipsForAtom(atom1)
1348	>	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1349	>	skip_it = .true.
1350	>	return
1351	>	endif
1352	>	end do
1353	>
1354	>	return
1355	>	end function skipThisPair
1356	>
1357	>	function FF_UsesDirectionalAtoms() result(doesit)
1358	>	logical :: doesit
1359	>	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
1362	>	end function FF_UsesDirectionalAtoms
1363	>
1364	>	function FF_RequiresPrepairCalc() result(doesit)
1365	>	logical :: doesit
1366	>	doesit = FF_uses_EAM
1367	>	end function FF_RequiresPrepairCalc
1368	>
1369	>	function FF_RequiresPostpairCalc() result(doesit)
1370	>	logical :: doesit
1371	>	doesit = FF_uses_RF
1372	>	end function FF_RequiresPostpairCalc
1373	>
1374		#ifdef PROFILE
1375	<	function getforcetime() result(totalforcetime)
1376	<	real(kind=dp) :: totalforcetime
1377	<	totalforcetime = forcetime
1378	<	end function getforcetime
1375	>	function getforcetime() result(totalforcetime)
1376	>	real(kind=dp) :: totalforcetime
1377	>	totalforcetime = forcetime
1378	>	end function getforcetime
1379		#endif
1158	–
1159	–	!! This cleans componets of force arrays belonging only to fortran
1380
1381	<	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
1381	>	!! This cleans componets of force arrays belonging only to fortran
1382
1383	<	!! Interfaces for C programs to module....
1383	>	subroutine add_stress_tensor(dpair, fpair)
1384
1385	<	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
1385	>	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1386
1387	<	integer, intent(out) :: thisStat
1388	<	call init_FF(LJMIXPOLICY, use_RF_c, thisStat)
1387	>	! because the d vector is the rj - ri vector, and
1388	>	! because fx, fy, fz are the force on atom i, we need a
1389	>	! negative sign here:
1390
1391	<	end subroutine initFortranFF
1391	>	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1392	>	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1393	>	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1394	>	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1395	>	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1396	>	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1397	>	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1398	>	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1399	>	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1400
1401	<	subroutine doForceloop(q, q_group, A, u_l, f, t, tau, pot, &
1402	<	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
1401	>	virial_Temp = virial_Temp + &
1402	>	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1403
1404	<	!! Stress Tensor
1405	<	real( kind = dp), dimension(9) :: tau
1406	<	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
1404	>	end subroutine add_stress_tensor
1405	>
1406	>	end module doForces

Diff Legend

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