ViewVC Help

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

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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 1610 by gezelter, Wed Oct 20 04:19:55 2004 UTC vs.
Revision 2262 by chuckv, Sun Jul 3 20:53:43 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.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
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
154		integer :: i
155	<	logical :: thisProperty
156	<	real (kind=DP) :: thisDPproperty
157	<
158	<	status = 0
159	<
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	>
175	>	if (.not. associated(atypes)) then
176	>	call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
177	>	status = -1
178	>	return
179	>	endif
180	>
181		nAtypes = getSize(atypes)
182	<
182	>
183		if (nAtypes == 0) then
184		status = -1
185		return
186		end if
119	–
120	–	if (.not. allocated(PropertyMap)) then
121	–	allocate(PropertyMap(nAtypes))
122	–	endif
187
188	+	if (.not. allocated(InteractionMap)) then
189	+	allocate(InteractionMap(nAtypes,nAtypes))
190	+	endif
191	+
192		do i = 1, nAtypes
193	<	call getElementProperty(atypes, i, "is_LJ", thisProperty)
194	<	PropertyMap(i)%is_LJ = thisProperty
195	<
196	<	call getElementProperty(atypes, i, "is_Charge", thisProperty)
197	<	PropertyMap(i)%is_Charge = thisProperty
198	<
199	<	if (thisProperty) then
132	<	call getElementProperty(atypes, i, "charge", thisDPproperty)
133	<	PropertyMap(i)%charge = thisDPproperty
134	<	endif
193	>	call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
194	>	call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
195	>	call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
196	>	call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
197	>	call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
198	>	call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
199	>	call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
200
201	<	call getElementProperty(atypes, i, "is_DP", thisProperty)
137	<	PropertyMap(i)%is_DP = thisProperty
201	>	do j = i, nAtypes
202
203	<	if (thisProperty) then
204	<	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
141	<	PropertyMap(i)%dipole_moment = thisDPproperty
142	<	endif
203	>	iHash = 0
204	>	myRcut = 0.0_dp
205
206	<	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
207	<	PropertyMap(i)%is_Sticky = thisProperty
208	<	call getElementProperty(atypes, i, "is_GB", thisProperty)
209	<	PropertyMap(i)%is_GB = thisProperty
210	<	call getElementProperty(atypes, i, "is_EAM", thisProperty)
211	<	PropertyMap(i)%is_EAM = thisProperty
206	>	call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
207	>	call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
208	>	call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
209	>	call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
210	>	call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
211	>	call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
212	>	call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
213	>
214	>	if (i_is_LJ .and. j_is_LJ) then
215	>	iHash = ior(iHash, LJ_PAIR)
216	>	endif
217	>
218	>	if (i_is_Elect .and. j_is_Elect) then
219	>	iHash = ior(iHash, ELECTROSTATIC_PAIR)
220	>	endif
221	>
222	>	if (i_is_Sticky .and. j_is_Sticky) then
223	>	iHash = ior(iHash, STICKY_PAIR)
224	>	endif
225	>
226	>	if (i_is_StickyP .and. j_is_StickyP) then
227	>	iHash = ior(iHash, STICKYPOWER_PAIR)
228	>	endif
229	>
230	>	if (i_is_EAM .and. j_is_EAM) then
231	>	iHash = ior(iHash, EAM_PAIR)
232	>	endif
233	>
234	>	if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
235	>	if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
236	>	if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
237	>
238	>	if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
239	>	if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
240	>	if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
241	>
242	>
243	>	InteractionMap(i,j)%InteractionHash = iHash
244	>	InteractionMap(j,i)%InteractionHash = iHash
245	>
246	>	end do
247	>
248		end do
249	+	end subroutine createInteractionMap
250
251	<	havePropertyMap = .true.
251	>	! 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.
252	>	subroutine createRcuts(defaultRList)
253	>	real(kind=dp), intent(in), optional :: defaultRList
254	>	integer :: iMap
255	>	integer :: map_i,map_j
256	>	real(kind=dp) :: thisRCut = 0.0_dp
257	>	real(kind=dp) :: actualCutoff = 0.0_dp
258	>	integer :: nAtypes
259
260	<	end subroutine createPropertyMap
260	>	if(.not. allocated(InteractionMap)) return
261
262	+	nAtypes = getSize(atypes)
263	+	! If we pass a default rcut, set all atypes to that cutoff distance
264	+	if(present(defaultRList)) then
265	+	InteractionMap(:,:)%rList = defaultRList
266	+	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
267	+	haveRlist = .true.
268	+	return
269	+	end if
270	+
271	+	do map_i = 1,nAtypes
272	+	do map_j = map_i,nAtypes
273	+	iMap = InteractionMap(map_i, map_j)%InteractionHash
274	+
275	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
276	+	!            thisRCut = getLJCutOff(map_i,map_j)
277	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
278	+	endif
279	+
280	+	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
281	+	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
282	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
283	+	endif
284	+
285	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
286	+	!             thisRCut = getStickyCutOff(map_i,map_j)
287	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
288	+	endif
289	+
290	+	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
291	+	!              thisRCut = getStickyPowerCutOff(map_i,map_j)
292	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
293	+	endif
294	+
295	+	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
296	+	!              thisRCut = getGayberneCutOff(map_i,map_j)
297	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
298	+	endif
299	+
300	+	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
301	+	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
302	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
303	+	endif
304	+
305	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
306	+	!              thisRCut = getEAMCutOff(map_i,map_j)
307	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
308	+	endif
309	+
310	+	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
311	+	!              thisRCut = getShapeCutOff(map_i,map_j)
312	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
313	+	endif
314	+
315	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
316	+	!              thisRCut = getShapeLJCutOff(map_i,map_j)
317	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
318	+	endif
319	+	InteractionMap(map_i, map_j)%rList = actualCutoff
320	+	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
321	+	end do
322	+	end do
323	+	haveRlist = .true.
324	+	end subroutine createRcuts
325	+
326	+
327	+	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
328	+	!!$  subroutine setRlistDF( this_rlist )
329	+	!!$
330	+	!!$   real(kind=dp) :: this_rlist
331	+	!!$
332	+	!!$    rlist = this_rlist
333	+	!!$    rlistsq = rlist * rlist
334	+	!!$
335	+	!!$    haveRlist = .true.
336	+	!!$
337	+	!!$  end subroutine setRlistDF
338	+
339	+
340		subroutine setSimVariables()
341	<	SIM_uses_LJ = SimUsesLJ()
342	<	SIM_uses_sticky = SimUsesSticky()
343	<	SIM_uses_charges = SimUsesCharges()
344	<	SIM_uses_dipoles = SimUsesDipoles()
345	<	SIM_uses_RF = SimUsesRF()
346	<	SIM_uses_GB = SimUsesGB()
341	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
342	>	SIM_uses_LennardJones = SimUsesLennardJones()
343	>	SIM_uses_Electrostatics = SimUsesElectrostatics()
344	>	SIM_uses_Charges = SimUsesCharges()
345	>	SIM_uses_Dipoles = SimUsesDipoles()
346	>	SIM_uses_Sticky = SimUsesSticky()
347	>	SIM_uses_StickyPower = SimUsesStickyPower()
348	>	SIM_uses_GayBerne = SimUsesGayBerne()
349		SIM_uses_EAM = SimUsesEAM()
350	+	SIM_uses_Shapes = SimUsesShapes()
351	+	SIM_uses_FLARB = SimUsesFLARB()
352	+	SIM_uses_RF = SimUsesRF()
353		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
354		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
166	–	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
355		SIM_uses_PBC = SimUsesPBC()
168	–	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
356
357		haveSIMvariables = .true.
358
#	Line 178 | Line 365 | contains
365		integer :: myStatus
366
367		error = 0
181	–
182	–	if (.not. havePropertyMap) then
368
369	+	if (.not. haveInteractionMap) then
370	+
371		myStatus = 0
372
373	<	call createPropertyMap(myStatus)
373	>	call createInteractionMap(myStatus)
374
375		if (myStatus .ne. 0) then
376	<	write(default_error, *) 'createPropertyMap failed in doForces!'
376	>	write(default_error, *) 'createInteractionMap failed in doForces!'
377		error = -1
378		return
379		endif
#	Line 202 | Line 389 | contains
389		return
390		endif
391
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	–
392		if (.not. haveNeighborList) then
393		write(default_error, *) 'neighbor list has not been initialized in doForces!'
394		error = -1
#	Line 231 | Line 410 | contains
410		#endif
411		return
412		end subroutine doReadyCheck
234	–
413
236	–	subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
414
415	<	integer, intent(in) :: LJMIXPOLICY
415	>	subroutine init_FF(use_RF_c, thisStat)
416	>
417		logical, intent(in) :: use_RF_c
418
419		integer, intent(out) :: thisStat
#	Line 248 | Line 426 | contains
426
427		!! Fortran's version of a cast:
428		FF_uses_RF = use_RF_c
429	<
429	>
430		!! init_FF is called *after* all of the atom types have been
431		!! defined in atype_module using the new_atype subroutine.
432		!!
433		!! this will scan through the known atypes and figure out what
434		!! interactions are used by the force field.
435	<
436	<	FF_uses_LJ = .false.
437	<	FF_uses_sticky = .false.
438	<	FF_uses_charges = .false.
439	<	FF_uses_dipoles = .false.
440	<	FF_uses_GB = .false.
435	>
436	>	FF_uses_DirectionalAtoms = .false.
437	>	FF_uses_LennardJones = .false.
438	>	FF_uses_Electrostatics = .false.
439	>	FF_uses_Charges = .false.
440	>	FF_uses_Dipoles = .false.
441	>	FF_uses_Sticky = .false.
442	>	FF_uses_StickyPower = .false.
443	>	FF_uses_GayBerne = .false.
444		FF_uses_EAM = .false.
445	<
446	<	call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
266	<	if (nMatches .gt. 0) FF_uses_LJ = .true.
445	>	FF_uses_Shapes = .false.
446	>	FF_uses_FLARB = .false.
447
448	<	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
449	<	if (nMatches .gt. 0) FF_uses_charges = .true.
448	>	call getMatchingElementList(atypes, "is_Directional", .true., &
449	>	nMatches, MatchList)
450	>	if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
451
452	<	call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
453	<	if (nMatches .gt. 0) FF_uses_dipoles = .true.
454	<
452	>	call getMatchingElementList(atypes, "is_LennardJones", .true., &
453	>	nMatches, MatchList)
454	>	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
455	>
456	>	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
457	>	nMatches, MatchList)
458	>	if (nMatches .gt. 0) then
459	>	FF_uses_Electrostatics = .true.
460	>	endif
461	>
462	>	call getMatchingElementList(atypes, "is_Charge", .true., &
463	>	nMatches, MatchList)
464	>	if (nMatches .gt. 0) then
465	>	FF_uses_Charges = .true.
466	>	FF_uses_Electrostatics = .true.
467	>	endif
468	>
469	>	call getMatchingElementList(atypes, "is_Dipole", .true., &
470	>	nMatches, MatchList)
471	>	if (nMatches .gt. 0) then
472	>	FF_uses_Dipoles = .true.
473	>	FF_uses_Electrostatics = .true.
474	>	FF_uses_DirectionalAtoms = .true.
475	>	endif
476	>
477	>	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
478	>	nMatches, MatchList)
479	>	if (nMatches .gt. 0) then
480	>	FF_uses_Quadrupoles = .true.
481	>	FF_uses_Electrostatics = .true.
482	>	FF_uses_DirectionalAtoms = .true.
483	>	endif
484	>
485		call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
486		MatchList)
487	<	if (nMatches .gt. 0) FF_uses_Sticky = .true.
487	>	if (nMatches .gt. 0) then
488	>	FF_uses_Sticky = .true.
489	>	FF_uses_DirectionalAtoms = .true.
490	>	endif
491	>
492	>	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
493	>	MatchList)
494	>	if (nMatches .gt. 0) then
495	>	FF_uses_StickyPower = .true.
496	>	FF_uses_DirectionalAtoms = .true.
497	>	endif
498
499	<	call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
500	<	if (nMatches .gt. 0) FF_uses_GB = .true.
501	<
499	>	call getMatchingElementList(atypes, "is_GayBerne", .true., &
500	>	nMatches, MatchList)
501	>	if (nMatches .gt. 0) then
502	>	FF_uses_GayBerne = .true.
503	>	FF_uses_DirectionalAtoms = .true.
504	>	endif
505	>
506		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
507		if (nMatches .gt. 0) FF_uses_EAM = .true.
508	<
508	>
509	>	call getMatchingElementList(atypes, "is_Shape", .true., &
510	>	nMatches, MatchList)
511	>	if (nMatches .gt. 0) then
512	>	FF_uses_Shapes = .true.
513	>	FF_uses_DirectionalAtoms = .true.
514	>	endif
515	>
516	>	call getMatchingElementList(atypes, "is_FLARB", .true., &
517	>	nMatches, MatchList)
518	>	if (nMatches .gt. 0) FF_uses_FLARB = .true.
519	>
520		!! Assume sanity (for the sake of argument)
521		haveSaneForceField = .true.
522
523		!! check to make sure the FF_uses_RF setting makes sense
524	<
524	>
525		if (FF_uses_dipoles) then
526		if (FF_uses_RF) then
527		dielect = getDielect()
#	Line 298 | Line 534 | contains
534		haveSaneForceField = .false.
535		return
536		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.
537		endif
538
539	<	if (FF_uses_sticky) then
540	<	call check_sticky_FF(my_status)
541	<	if (my_status /= 0) then
542	<	thisStat = -1
543	<	haveSaneForceField = .false.
544	<	return
545	<	end if
546	<	endif
539	>	!sticky module does not contain check_sticky_FF anymore
540	>	!if (FF_uses_sticky) then
541	>	!   call check_sticky_FF(my_status)
542	>	!   if (my_status /= 0) then
543	>	!      thisStat = -1
544	>	!      haveSaneForceField = .false.
545	>	!      return
546	>	!   end if
547	>	!endif
548
333	–
549		if (FF_uses_EAM) then
550	<	call init_EAM_FF(my_status)
550	>	call init_EAM_FF(my_status)
551		if (my_status /= 0) then
552		write(default_error, *) "init_EAM_FF returned a bad status"
553		thisStat = -1
#	Line 341 | Line 556 | contains
556		end if
557		endif
558
559	<	if (FF_uses_GB) then
559	>	if (FF_uses_GayBerne) then
560		call check_gb_pair_FF(my_status)
561		if (my_status .ne. 0) then
562		thisStat = -1
#	Line 350 | Line 565 | contains
565		endif
566		endif
567
568	<	if (FF_uses_GB .and. FF_uses_LJ) then
568	>	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569		endif
570	+
571		if (.not. haveNeighborList) then
572		!! Create neighbor lists
573		call expandNeighborList(nLocal, my_status)
#	Line 363 | Line 579 | contains
579		haveNeighborList = .true.
580		endif
581
366	–
367	–
582		end subroutine init_FF
369	–
583
584	+
585		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
586		!------------------------------------------------------------->
587	<	subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
587	>	subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
588		do_pot_c, do_stress_c, error)
589		!! Position array provided by C, dimensioned by getNlocal
590		real ( kind = dp ), dimension(3, nLocal) :: q
#	Line 379 | Line 593 | contains
593		!! Rotation Matrix for each long range particle in simulation.
594		real( kind = dp), dimension(9, nLocal) :: A
595		!! Unit vectors for dipoles (lab frame)
596	<	real( kind = dp ), dimension(3,nLocal) :: u_l
596	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
597		!! Force array provided by C, dimensioned by getNlocal
598		real ( kind = dp ), dimension(3,nLocal) :: f
599		!! Torsion array provided by C, dimensioned by getNlocal
#	Line 417 | Line 631 | contains
631		integer :: localError
632		integer :: propPack_i, propPack_j
633		integer :: loopStart, loopEnd, loop
634	<
634	>	integer :: iMap
635		real(kind=dp) :: listSkin = 1.0
636	<
636	>
637		!! initialize local variables
638	<
638	>
639		#ifdef IS_MPI
640		pot_local = 0.0_dp
641		nAtomsInRow   = getNatomsInRow(plan_atom_row)
#	Line 431 | Line 645 | contains
645		#else
646		natoms = nlocal
647		#endif
648	<
648	>
649		call doReadyCheck(localError)
650		if ( localError .ne. 0 ) then
651		call handleError("do_force_loop", "Not Initialized")
#	Line 439 | Line 653 | contains
653		return
654		end if
655		call zero_work_arrays()
656	<
656	>
657		do_pot = do_pot_c
658		do_stress = do_stress_c
659	<
659	>
660		! Gather all information needed by all force loops:
661	<
661	>
662		#ifdef IS_MPI
663	<
663	>
664		call gather(q, q_Row, plan_atom_row_3d)
665		call gather(q, q_Col, plan_atom_col_3d)
666
667		call gather(q_group, q_group_Row, plan_group_row_3d)
668		call gather(q_group, q_group_Col, plan_group_col_3d)
669	<
670	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
671	<	call gather(u_l, u_l_Row, plan_atom_row_3d)
672	<	call gather(u_l, u_l_Col, plan_atom_col_3d)
673	<
669	>
670	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
671	>	call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
672	>	call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
673	>
674		call gather(A, A_Row, plan_atom_row_rotation)
675		call gather(A, A_Col, plan_atom_col_rotation)
676		endif
677	<
677	>
678		#endif
679	<
679	>
680		!! Begin force loop timing:
681		#ifdef PROFILE
682		call cpu_time(forceTimeInitial)
683		nloops = nloops + 1
684		#endif
685	<
685	>
686		loopEnd = PAIR_LOOP
687		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
688		loopStart = PREPAIR_LOOP
#	Line 483 | Line 697 | contains
697		if (loop .eq. loopStart) then
698		#ifdef IS_MPI
699		call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
700	<	update_nlist)
700	>	update_nlist)
701		#else
702		call checkNeighborList(nGroups, q_group, listSkin, &
703	<	update_nlist)
703	>	update_nlist)
704		#endif
705		endif
706	<
706	>
707		if (update_nlist) then
708		!! save current configuration and construct neighbor list
709		#ifdef IS_MPI
#	Line 500 | Line 714 | contains
714		neighborListSize = size(list)
715		nlist = 0
716		endif
717	<
717	>
718		istart = 1
719		#ifdef IS_MPI
720		iend = nGroupsInRow
#	Line 510 | Line 724 | contains
724		outer: do i = istart, iend
725
726		if (update_nlist) point(i) = nlist + 1
727	<
727	>
728		n_in_i = groupStartRow(i+1) - groupStartRow(i)
729	<
729	>
730		if (update_nlist) then
731		#ifdef IS_MPI
732		jstart = 1
#	Line 527 | Line 741 | contains
741		! make sure group i has neighbors
742		if (jstart .gt. jend) cycle outer
743		endif
744	<
744	>
745		do jnab = jstart, jend
746		if (update_nlist) then
747		j = jnab
#	Line 543 | Line 757 | contains
757		q_group(:,j), d_grp, rgrpsq)
758		#endif
759
760	<	if (rgrpsq < rlistsq) then
760	>	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
761		if (update_nlist) then
762		nlist = nlist + 1
763	<
763	>
764		if (nlist > neighborListSize) then
765		#ifdef IS_MPI
766		call expandNeighborList(nGroupsInRow, listerror)
#	Line 560 | Line 774 | contains
774		end if
775		neighborListSize = size(list)
776		endif
777	<
777	>
778		list(nlist) = j
779		endif
780	<
780	>
781		if (loop .eq. PAIR_LOOP) then
782		vij = 0.0d0
783		fij(1:3) = 0.0d0
784		endif
785	<
785	>
786		call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
787		in_switching_region)
788	<
788	>
789		n_in_j = groupStartCol(j+1) - groupStartCol(j)
790	<
790	>
791		do ia = groupStartRow(i), groupStartRow(i+1)-1
792	<
792	>
793		atom1 = groupListRow(ia)
794	<
794	>
795		inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
796	<
796	>
797		atom2 = groupListCol(jb)
798	<
798	>
799		if (skipThisPair(atom1, atom2)) cycle inner
800
801		if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
#	Line 601 | Line 815 | contains
815		#ifdef IS_MPI
816		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
817		rgrpsq, d_grp, do_pot, do_stress, &
818	<	u_l, A, f, t, pot_local)
818	>	eFrame, A, f, t, pot_local)
819		#else
820		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
821		rgrpsq, d_grp, do_pot, do_stress, &
822	<	u_l, A, f, t, pot)
822	>	eFrame, A, f, t, pot)
823		#endif
824		else
825		#ifdef IS_MPI
826		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
827		do_pot, &
828	<	u_l, A, f, t, pot_local, vpair, fpair)
828	>	eFrame, A, f, t, pot_local, vpair, fpair)
829		#else
830		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
831		do_pot,  &
832	<	u_l, A, f, t, pot, vpair, fpair)
832	>	eFrame, A, f, t, pot, vpair, fpair)
833		#endif
834
835		vij = vij + vpair
#	Line 623 | Line 837 | contains
837		endif
838		enddo inner
839		enddo
840	<
840	>
841		if (loop .eq. PAIR_LOOP) then
842		if (in_switching_region) then
843		swderiv = vij*dswdr/rgrp
844		fij(1) = fij(1) + swderiv*d_grp(1)
845		fij(2) = fij(2) + swderiv*d_grp(2)
846		fij(3) = fij(3) + swderiv*d_grp(3)
847	<
847	>
848		do ia=groupStartRow(i), groupStartRow(i+1)-1
849		atom1=groupListRow(ia)
850		mf = mfactRow(atom1)
#	Line 644 | Line 858 | contains
858		f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
859		#endif
860		enddo
861	<
861	>
862		do jb=groupStartCol(j), groupStartCol(j+1)-1
863		atom2=groupListCol(jb)
864		mf = mfactCol(atom2)
#	Line 659 | Line 873 | contains
873		#endif
874		enddo
875		endif
876	<
876	>
877		if (do_stress) call add_stress_tensor(d_grp, fij)
878		endif
879		end if
880		enddo
881		enddo outer
882	<
882	>
883		if (update_nlist) then
884		#ifdef IS_MPI
885		point(nGroupsInRow + 1) = nlist + 1
#	Line 679 | Line 893 | contains
893		update_nlist = .false.
894		endif
895		endif
896	<
896	>
897		if (loop .eq. PREPAIR_LOOP) then
898		call do_preforce(nlocal, pot)
899		endif
900	<
900	>
901		enddo
902	<
902	>
903		!! Do timing
904		#ifdef PROFILE
905		call cpu_time(forceTimeFinal)
906		forceTime = forceTime + forceTimeFinal - forceTimeInitial
907		#endif
908	<
908	>
909		#ifdef IS_MPI
910		!!distribute forces
911	<
911	>
912		f_temp = 0.0_dp
913		call scatter(f_Row,f_temp,plan_atom_row_3d)
914		do i = 1,nlocal
915		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
916		end do
917	<
917	>
918		f_temp = 0.0_dp
919		call scatter(f_Col,f_temp,plan_atom_col_3d)
920		do i = 1,nlocal
921		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
922		end do
923	<
924	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
923	>
924	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
925		t_temp = 0.0_dp
926		call scatter(t_Row,t_temp,plan_atom_row_3d)
927		do i = 1,nlocal
#	Line 715 | Line 929 | contains
929		end do
930		t_temp = 0.0_dp
931		call scatter(t_Col,t_temp,plan_atom_col_3d)
932	<
932	>
933		do i = 1,nlocal
934		t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
935		end do
936		endif
937	<
937	>
938		if (do_pot) then
939		! scatter/gather pot_row into the members of my column
940		call scatter(pot_Row, pot_Temp, plan_atom_row)
941	<
941	>
942		! scatter/gather pot_local into all other procs
943		! add resultant to get total pot
944		do i = 1, nlocal
945		pot_local = pot_local + pot_Temp(i)
946		enddo
947	<
947	>
948		pot_Temp = 0.0_DP
949	<
949	>
950		call scatter(pot_Col, pot_Temp, plan_atom_col)
951		do i = 1, nlocal
952		pot_local = pot_local + pot_Temp(i)
953		enddo
954	<
954	>
955		endif
956		#endif
957	<
957	>
958		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
959	<
959	>
960		if (FF_uses_RF .and. SIM_uses_RF) then
961	<
961	>
962		#ifdef IS_MPI
963		call scatter(rf_Row,rf,plan_atom_row_3d)
964		call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
#	Line 752 | Line 966 | contains
966		rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
967		end do
968		#endif
969	<
969	>
970		do i = 1, nLocal
971	<
971	>
972		rfpot = 0.0_DP
973		#ifdef IS_MPI
974		me_i = atid_row(i)
975		#else
976		me_i = atid(i)
977		#endif
978	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
979
980	<	if (PropertyMap(me_i)%is_DP) then
981	<
982	<	mu_i = PropertyMap(me_i)%dipole_moment
983	<
980	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
981	>
982	>	mu_i = getDipoleMoment(me_i)
983	>
984		!! The reaction field needs to include a self contribution
985		!! to the field:
986	<	call accumulate_self_rf(i, mu_i, u_l)
986	>	call accumulate_self_rf(i, mu_i, eFrame)
987		!! Get the reaction field contribution to the
988		!! potential and torques:
989	<	call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
989	>	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
990		#ifdef IS_MPI
991		pot_local = pot_local + rfpot
992		#else
993		pot = pot + rfpot
994	<
994	>
995		#endif
996	<	endif
996	>	endif
997		enddo
998		endif
999		endif
1000	<
1001	<
1000	>
1001	>
1002		#ifdef IS_MPI
1003	<
1003	>
1004		if (do_pot) then
1005		pot = pot + pot_local
1006		!! we assume the c code will do the allreduce to get the total potential
1007		!! we could do it right here if we needed to...
1008		endif
1009	<
1009	>
1010		if (do_stress) then
1011		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1012		mpi_comm_world,mpi_err)
1013		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1014		mpi_comm_world,mpi_err)
1015		endif
1016	<
1016	>
1017		#else
1018	<
1018	>
1019		if (do_stress) then
1020		tau = tau_Temp
1021		virial = virial_Temp
1022		endif
1023	<
1023	>
1024		#endif
1025	<
1025	>
1026		end subroutine do_force_loop
1027	<
1027	>
1028		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1029	<	u_l, A, f, t, pot, vpair, fpair)
1029	>	eFrame, A, f, t, pot, vpair, fpair)
1030
1031		real( kind = dp ) :: pot, vpair, sw
1032		real( kind = dp ), dimension(3) :: fpair
1033		real( kind = dp ), dimension(nLocal)   :: mfact
1034	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1034	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1035		real( kind = dp ), dimension(9,nLocal) :: A
1036		real( kind = dp ), dimension(3,nLocal) :: f
1037		real( kind = dp ), dimension(3,nLocal) :: t
#	Line 826 | Line 1041 | contains
1041		real ( kind = dp ), intent(inout) :: rijsq
1042		real ( kind = dp )                :: r
1043		real ( kind = dp ), intent(inout) :: d(3)
1044	+	real ( kind = dp ) :: ebalance
1045		integer :: me_i, me_j
1046
1047	+	integer :: iMap
1048	+
1049		r = sqrt(rijsq)
1050		vpair = 0.0d0
1051		fpair(1:3) = 0.0d0
#	Line 839 | Line 1057 | contains
1057		me_i = atid(i)
1058		me_j = atid(j)
1059		#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
1060
1061	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
1062	+
1063	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1064	+	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1065		endif
1066
1067	<	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1067	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1068	>	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1069	>	pot, eFrame, f, t, do_pot)
1070
1071	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1072	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
1073	<	do_pot)
1071	>	if (FF_uses_RF .and. SIM_uses_RF) then
1072	>
1073	>	! CHECK ME (RF needs to know about all electrostatic types)
1074	>	call accumulate_rf(i, j, r, eFrame, sw)
1075	>	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1076		endif
1077
1078		endif
1079
1080	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1081	+	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1082	+	pot, A, f, t, do_pot)
1083	+	endif
1084
1085	<	if (FF_uses_GB .and. SIM_uses_GB) then
1086	<
1087	<	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1088	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
886	<	do_pot)
887	<	endif
1085	>	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1086	>	call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087	>	pot, A, f, t, do_pot)
1088	>	endif
1089
1090	+	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1091	+	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1092	+	pot, A, f, t, do_pot)
1093		endif
1094	<
1095	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1096	<
1097	<	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	<
1094	>
1095	>	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1096	>	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1097	>	!           pot, A, f, t, do_pot)
1098		endif
1099	+
1100	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1101	+	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1102	+	do_pot)
1103	+	endif
1104	+
1105	+	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1106	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1107	+	pot, A, f, t, do_pot)
1108	+	endif
1109	+
1110	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1111	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1112	+	pot, A, f, t, do_pot)
1113	+	endif
1114
1115		end subroutine do_pair
1116
1117		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1118	<	do_pot, do_stress, u_l, A, f, t, pot)
1118	>	do_pot, do_stress, eFrame, A, f, t, pot)
1119
1120	<	real( kind = dp ) :: pot, sw
1121	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1122	<	real (kind=dp), dimension(9,nLocal) :: A
1123	<	real (kind=dp), dimension(3,nLocal) :: f
1124	<	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	<
1120	>	real( kind = dp ) :: pot, sw
1121	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1122	>	real (kind=dp), dimension(9,nLocal) :: A
1123	>	real (kind=dp), dimension(3,nLocal) :: f
1124	>	real (kind=dp), dimension(3,nLocal) :: t
1125
1126	<	r = sqrt(rijsq)
1127	<	if (SIM_uses_molecular_cutoffs) then
1128	<	rc = sqrt(rcijsq)
1129	<	else
1130	<	rc = r
927	<	endif
928	<
1126	>	logical, intent(inout) :: do_pot, do_stress
1127	>	integer, intent(in) :: i, j
1128	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1129	>	real ( kind = dp )                :: r, rc
1130	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1131
1132	+	integer :: me_i, me_j, iMap
1133	+
1134		#ifdef IS_MPI
1135	<	me_i = atid_row(i)
1136	<	me_j = atid_col(j)
1135	>	me_i = atid_row(i)
1136	>	me_j = atid_col(j)
1137		#else
1138	<	me_i = atid(i)
1139	<	me_j = atid(j)
1138	>	me_i = atid(i)
1139	>	me_j = atid(j)
1140		#endif
1141	<
1142	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1143	<
1144	<	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1145	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1146	<
1147	<	endif
1148	<
1149	<	end subroutine do_prepair
1150	<
1151	<
1152	<	subroutine do_preforce(nlocal,pot)
1153	<	integer :: nlocal
1154	<	real( kind = dp ) :: pot
1155	<
1156	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1157	<	call calc_EAM_preforce_Frho(nlocal,pot)
1158	<	endif
1159	<
1160	<
1161	<	end subroutine do_preforce
1162	<
1163	<
1164	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1165	<
1166	<	real (kind = dp), dimension(3) :: q_i
1167	<	real (kind = dp), dimension(3) :: q_j
1168	<	real ( kind = dp ), intent(out) :: r_sq
1169	<	real( kind = dp ) :: d(3), scaled(3)
1170	<	integer i
1171	<
1172	<	d(1:3) = q_j(1:3) - q_i(1:3)
1173	<
1174	<	! Wrap back into periodic box if necessary
1175	<	if ( SIM_uses_PBC ) then
1176	<
1177	<	if( .not.boxIsOrthorhombic ) then
1178	<	! calc the scaled coordinates.
1179	<
1180	<	scaled = matmul(HmatInv, d)
1181	<
1182	<	! wrap the scaled coordinates
1183	<
1184	<	scaled = scaled  - anint(scaled)
1185	<
1186	<
1187	<	! calc the wrapped real coordinates from the wrapped scaled
1188	<	! coordinates
1189	<
1190	<	d = matmul(Hmat,scaled)
1191	<
1192	<	else
1193	<	! calc the scaled coordinates.
1194	<
1195	<	do i = 1, 3
1196	<	scaled(i) = d(i) * HmatInv(i,i)
1197	<
1198	<	! wrap the scaled coordinates
1199	<
1200	<	scaled(i) = scaled(i) - anint(scaled(i))
1201	<
1202	<	! calc the wrapped real coordinates from the wrapped scaled
1203	<	! coordinates
1204	<
1205	<	d(i) = scaled(i)*Hmat(i,i)
1206	<	enddo
1207	<	endif
1208	<
1209	<	endif
1210	<
1211	<	r_sq = dot_product(d,d)
1212	<
1213	<	end subroutine get_interatomic_vector
1214	<
1215	<	subroutine zero_work_arrays()
1012	<
1141	>
1142	>	iMap = InteractionMap(me_i, me_j)%InteractionHash
1143	>
1144	>	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1145	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1146	>	endif
1147	>
1148	>	end subroutine do_prepair
1149	>
1150	>
1151	>	subroutine do_preforce(nlocal,pot)
1152	>	integer :: nlocal
1153	>	real( kind = dp ) :: pot
1154	>
1155	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1156	>	call calc_EAM_preforce_Frho(nlocal,pot)
1157	>	endif
1158	>
1159	>
1160	>	end subroutine do_preforce
1161	>
1162	>
1163	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1164	>
1165	>	real (kind = dp), dimension(3) :: q_i
1166	>	real (kind = dp), dimension(3) :: q_j
1167	>	real ( kind = dp ), intent(out) :: r_sq
1168	>	real( kind = dp ) :: d(3), scaled(3)
1169	>	integer i
1170	>
1171	>	d(1:3) = q_j(1:3) - q_i(1:3)
1172	>
1173	>	! Wrap back into periodic box if necessary
1174	>	if ( SIM_uses_PBC ) then
1175	>
1176	>	if( .not.boxIsOrthorhombic ) then
1177	>	! calc the scaled coordinates.
1178	>
1179	>	scaled = matmul(HmatInv, d)
1180	>
1181	>	! wrap the scaled coordinates
1182	>
1183	>	scaled = scaled  - anint(scaled)
1184	>
1185	>
1186	>	! calc the wrapped real coordinates from the wrapped scaled
1187	>	! coordinates
1188	>
1189	>	d = matmul(Hmat,scaled)
1190	>
1191	>	else
1192	>	! calc the scaled coordinates.
1193	>
1194	>	do i = 1, 3
1195	>	scaled(i) = d(i) * HmatInv(i,i)
1196	>
1197	>	! wrap the scaled coordinates
1198	>
1199	>	scaled(i) = scaled(i) - anint(scaled(i))
1200	>
1201	>	! calc the wrapped real coordinates from the wrapped scaled
1202	>	! coordinates
1203	>
1204	>	d(i) = scaled(i)*Hmat(i,i)
1205	>	enddo
1206	>	endif
1207	>
1208	>	endif
1209	>
1210	>	r_sq = dot_product(d,d)
1211	>
1212	>	end subroutine get_interatomic_vector
1213	>
1214	>	subroutine zero_work_arrays()
1215	>
1216		#ifdef IS_MPI
1014	–
1015	–	q_Row = 0.0_dp
1016	–	q_Col = 0.0_dp
1217
1218	<	q_group_Row = 0.0_dp
1219	<	q_group_Col = 0.0_dp
1220	<
1221	<	u_l_Row = 0.0_dp
1222	<	u_l_Col = 0.0_dp
1223	<
1224	<	A_Row = 0.0_dp
1225	<	A_Col = 0.0_dp
1226	<
1227	<	f_Row = 0.0_dp
1228	<	f_Col = 0.0_dp
1229	<	f_Temp = 0.0_dp
1230	<
1231	<	t_Row = 0.0_dp
1232	<	t_Col = 0.0_dp
1233	<	t_Temp = 0.0_dp
1234	<
1235	<	pot_Row = 0.0_dp
1236	<	pot_Col = 0.0_dp
1237	<	pot_Temp = 0.0_dp
1238	<
1239	<	rf_Row = 0.0_dp
1240	<	rf_Col = 0.0_dp
1241	<	rf_Temp = 0.0_dp
1242	<
1218	>	q_Row = 0.0_dp
1219	>	q_Col = 0.0_dp
1220	>
1221	>	q_group_Row = 0.0_dp
1222	>	q_group_Col = 0.0_dp
1223	>
1224	>	eFrame_Row = 0.0_dp
1225	>	eFrame_Col = 0.0_dp
1226	>
1227	>	A_Row = 0.0_dp
1228	>	A_Col = 0.0_dp
1229	>
1230	>	f_Row = 0.0_dp
1231	>	f_Col = 0.0_dp
1232	>	f_Temp = 0.0_dp
1233	>
1234	>	t_Row = 0.0_dp
1235	>	t_Col = 0.0_dp
1236	>	t_Temp = 0.0_dp
1237	>
1238	>	pot_Row = 0.0_dp
1239	>	pot_Col = 0.0_dp
1240	>	pot_Temp = 0.0_dp
1241	>
1242	>	rf_Row = 0.0_dp
1243	>	rf_Col = 0.0_dp
1244	>	rf_Temp = 0.0_dp
1245	>
1246		#endif
1247	<
1248	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1249	<	call clean_EAM()
1250	<	endif
1251	<
1252	<	rf = 0.0_dp
1253	<	tau_Temp = 0.0_dp
1254	<	virial_Temp = 0.0_dp
1255	<	end subroutine zero_work_arrays
1256	<
1257	<	function skipThisPair(atom1, atom2) result(skip_it)
1258	<	integer, intent(in) :: atom1
1259	<	integer, intent(in), optional :: atom2
1260	<	logical :: skip_it
1261	<	integer :: unique_id_1, unique_id_2
1262	<	integer :: me_i,me_j
1263	<	integer :: i
1264	<
1265	<	skip_it = .false.
1266	<
1267	<	!! there are a number of reasons to skip a pair or a particle
1268	<	!! mostly we do this to exclude atoms who are involved in short
1269	<	!! range interactions (bonds, bends, torsions), but we also need
1270	<	!! to exclude some overcounted interactions that result from
1271	<	!! the parallel decomposition
1272	<
1247	>
1248	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1249	>	call clean_EAM()
1250	>	endif
1251	>
1252	>	rf = 0.0_dp
1253	>	tau_Temp = 0.0_dp
1254	>	virial_Temp = 0.0_dp
1255	>	end subroutine zero_work_arrays
1256	>
1257	>	function skipThisPair(atom1, atom2) result(skip_it)
1258	>	integer, intent(in) :: atom1
1259	>	integer, intent(in), optional :: atom2
1260	>	logical :: skip_it
1261	>	integer :: unique_id_1, unique_id_2
1262	>	integer :: me_i,me_j
1263	>	integer :: i
1264	>
1265	>	skip_it = .false.
1266	>
1267	>	!! there are a number of reasons to skip a pair or a particle
1268	>	!! mostly we do this to exclude atoms who are involved in short
1269	>	!! range interactions (bonds, bends, torsions), but we also need
1270	>	!! to exclude some overcounted interactions that result from
1271	>	!! the parallel decomposition
1272	>
1273		#ifdef IS_MPI
1274	<	!! in MPI, we have to look up the unique IDs for each atom
1275	<	unique_id_1 = AtomRowToGlobal(atom1)
1274	>	!! in MPI, we have to look up the unique IDs for each atom
1275	>	unique_id_1 = AtomRowToGlobal(atom1)
1276		#else
1277	<	!! in the normal loop, the atom numbers are unique
1278	<	unique_id_1 = atom1
1277	>	!! in the normal loop, the atom numbers are unique
1278	>	unique_id_1 = atom1
1279		#endif
1280	<
1281	<	!! We were called with only one atom, so just check the global exclude
1282	<	!! list for this atom
1283	<	if (.not. present(atom2)) then
1284	<	do i = 1, nExcludes_global
1285	<	if (excludesGlobal(i) == unique_id_1) then
1286	<	skip_it = .true.
1287	<	return
1288	<	end if
1289	<	end do
1290	<	return
1291	<	end if
1292	<
1280	>
1281	>	!! We were called with only one atom, so just check the global exclude
1282	>	!! list for this atom
1283	>	if (.not. present(atom2)) then
1284	>	do i = 1, nExcludes_global
1285	>	if (excludesGlobal(i) == unique_id_1) then
1286	>	skip_it = .true.
1287	>	return
1288	>	end if
1289	>	end do
1290	>	return
1291	>	end if
1292	>
1293		#ifdef IS_MPI
1294	<	unique_id_2 = AtomColToGlobal(atom2)
1294	>	unique_id_2 = AtomColToGlobal(atom2)
1295		#else
1296	<	unique_id_2 = atom2
1296	>	unique_id_2 = atom2
1297		#endif
1298	<
1298	>
1299		#ifdef IS_MPI
1300	<	!! this situation should only arise in MPI simulations
1301	<	if (unique_id_1 == unique_id_2) then
1302	<	skip_it = .true.
1303	<	return
1304	<	end if
1305	<
1306	<	!! this prevents us from doing the pair on multiple processors
1307	<	if (unique_id_1 < unique_id_2) then
1308	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1309	<	skip_it = .true.
1310	<	return
1311	<	endif
1312	<	else
1313	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1314	<	skip_it = .true.
1315	<	return
1316	<	endif
1317	<	endif
1300	>	!! this situation should only arise in MPI simulations
1301	>	if (unique_id_1 == unique_id_2) then
1302	>	skip_it = .true.
1303	>	return
1304	>	end if
1305	>
1306	>	!! this prevents us from doing the pair on multiple processors
1307	>	if (unique_id_1 < unique_id_2) then
1308	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1309	>	skip_it = .true.
1310	>	return
1311	>	endif
1312	>	else
1313	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1314	>	skip_it = .true.
1315	>	return
1316	>	endif
1317	>	endif
1318		#endif
1319	<
1320	<	!! the rest of these situations can happen in all simulations:
1321	<	do i = 1, nExcludes_global
1322	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1323	<	(excludesGlobal(i) == unique_id_2)) then
1324	<	skip_it = .true.
1325	<	return
1326	<	endif
1327	<	enddo
1328	<
1329	<	do i = 1, nSkipsForAtom(atom1)
1330	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1331	<	skip_it = .true.
1332	<	return
1333	<	endif
1334	<	end do
1335	<
1336	<	return
1337	<	end function skipThisPair
1338	<
1339	<	function FF_UsesDirectionalAtoms() result(doesit)
1340	<	logical :: doesit
1341	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1342	<	FF_uses_GB .or. FF_uses_RF
1343	<	end function FF_UsesDirectionalAtoms
1344	<
1345	<	function FF_RequiresPrepairCalc() result(doesit)
1346	<	logical :: doesit
1347	<	doesit = FF_uses_EAM
1348	<	end function FF_RequiresPrepairCalc
1349	<
1350	<	function FF_RequiresPostpairCalc() result(doesit)
1351	<	logical :: doesit
1352	<	doesit = FF_uses_RF
1353	<	end function FF_RequiresPostpairCalc
1354	<
1319	>
1320	>	!! the rest of these situations can happen in all simulations:
1321	>	do i = 1, nExcludes_global
1322	>	if ((excludesGlobal(i) == unique_id_1) .or. &
1323	>	(excludesGlobal(i) == unique_id_2)) then
1324	>	skip_it = .true.
1325	>	return
1326	>	endif
1327	>	enddo
1328	>
1329	>	do i = 1, nSkipsForAtom(atom1)
1330	>	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1331	>	skip_it = .true.
1332	>	return
1333	>	endif
1334	>	end do
1335	>
1336	>	return
1337	>	end function skipThisPair
1338	>
1339	>	function FF_UsesDirectionalAtoms() result(doesit)
1340	>	logical :: doesit
1341	>	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1342	>	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1343	>	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1344	>	end function FF_UsesDirectionalAtoms
1345	>
1346	>	function FF_RequiresPrepairCalc() result(doesit)
1347	>	logical :: doesit
1348	>	doesit = FF_uses_EAM
1349	>	end function FF_RequiresPrepairCalc
1350	>
1351	>	function FF_RequiresPostpairCalc() result(doesit)
1352	>	logical :: doesit
1353	>	doesit = FF_uses_RF
1354	>	end function FF_RequiresPostpairCalc
1355	>
1356		#ifdef PROFILE
1357	<	function getforcetime() result(totalforcetime)
1358	<	real(kind=dp) :: totalforcetime
1359	<	totalforcetime = forcetime
1360	<	end function getforcetime
1357	>	function getforcetime() result(totalforcetime)
1358	>	real(kind=dp) :: totalforcetime
1359	>	totalforcetime = forcetime
1360	>	end function getforcetime
1361		#endif
1158	–
1159	–	!! This cleans componets of force arrays belonging only to fortran
1362
1363	<	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
1363	>	!! This cleans componets of force arrays belonging only to fortran
1364
1365	<	!! Interfaces for C programs to module....
1365	>	subroutine add_stress_tensor(dpair, fpair)
1366
1367	<	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
1367	>	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1368
1369	<	integer, intent(out) :: thisStat
1370	<	call init_FF(LJMIXPOLICY, use_RF_c, thisStat)
1369	>	! because the d vector is the rj - ri vector, and
1370	>	! because fx, fy, fz are the force on atom i, we need a
1371	>	! negative sign here:
1372
1373	<	end subroutine initFortranFF
1373	>	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1374	>	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1375	>	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1376	>	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1377	>	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1378	>	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1379	>	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1380	>	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1381	>	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1382
1383	<	subroutine doForceloop(q, q_group, A, u_l, f, t, tau, pot, &
1384	<	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
1383	>	virial_Temp = virial_Temp + &
1384	>	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1385
1386	<	!! Stress Tensor
1387	<	real( kind = dp), dimension(9) :: tau
1388	<	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
1386	>	end subroutine add_stress_tensor
1387	>
1388	>	end module doForces

Diff Legend

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