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 1628 by gezelter, Thu Oct 21 20:15:31 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.2 2004-10-21 20:15:22 gezelter Exp $, $Date: 2004-10-21 20:15:22 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
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 31 | Line 73 | module doForces
73
74		#define __FORTRAN90
75		#include "UseTheForce/fSwitchingFunction.h"
76	+	#include "UseTheForce/DarkSide/fInteractionMap.h"
77
78		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
79		INTEGER, PARAMETER:: PAIR_LOOP    = 2
#	Line 38 | Line 81 | module doForces
81		logical, save :: haveRlist = .false.
82		logical, save :: haveNeighborList = .false.
83		logical, save :: haveSIMvariables = .false.
41	–	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
59	–	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 73 | 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
84	<	real(kind=DP) :: dipole_moment = 0.0_DP
85	<	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
88	<
147	>
148		contains
149
150	<	subroutine setRlistDF( this_rlist )
151	<
93	<	real(kind=dp) :: this_rlist
94	<
95	<	rlist = this_rlist
96	<	rlistsq = rlist * rlist
97	<
98	<	haveRlist = .true.
99	<
100	<	end subroutine setRlistDF
101	<
102	<	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
117	–
118	–	if (.not. allocated(PropertyMap)) then
119	–	allocate(PropertyMap(nAtypes))
120	–	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)
127	<	PropertyMap(i)%is_Charge = thisProperty
128	<
129	<	if (thisProperty) then
130	<	call getElementProperty(atypes, i, "charge", thisDPproperty)
131	<	PropertyMap(i)%charge = thisDPproperty
132	<	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()
164	–	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
369		SIM_uses_PBC = SimUsesPBC()
166	–	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
370
371		haveSIMvariables = .true.
372
#	Line 176 | Line 379 | contains
379		integer :: myStatus
380
381		error = 0
179	–
180	–	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 221 | Line 424 | contains
424		#endif
425		return
426		end subroutine doReadyCheck
224	–
427
428	+
429		subroutine init_FF(use_RF_c, thisStat)
430
431		logical, intent(in) :: use_RF_c
#	Line 237 | 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)
461	<	if (nMatches .gt. 0) FF_uses_LJ = .true.
462	<
463	<	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
464	<	if (nMatches .gt. 0) FF_uses_charges = .true.
465	<
466	<	call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
467	<	if (nMatches .gt. 0) FF_uses_dipoles = .true.
468	<
459	>	FF_uses_Shapes = .false.
460	>	FF_uses_FLARB = .false.
461	>
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_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.
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_GB", .true., nMatches, MatchList)
514	<	if (nMatches .gt. 0) FF_uses_GB = .true.
515	<
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	<
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 287 | Line 548 | contains
548		haveSaneForceField = .false.
549		return
550		endif
290	–	endif
291	–
292	–	if (FF_uses_sticky) then
293	–	call check_sticky_FF(my_status)
294	–	if (my_status /= 0) then
295	–	thisStat = -1
296	–	haveSaneForceField = .false.
297	–	return
298	–	end if
551		endif
552
553	+	!sticky module does not contain check_sticky_FF anymore
554	+	!if (FF_uses_sticky) then
555	+	!   call check_sticky_FF(my_status)
556	+	!   if (my_status /= 0) then
557	+	!      thisStat = -1
558	+	!      haveSaneForceField = .false.
559	+	!      return
560	+	!   end if
561	+	!endif
562	+
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 308 | 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 317 | 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
#	Line 329 | Line 591 | contains
591		return
592		endif
593		haveNeighborList = .true.
594	<	endif
595	<
594	>	endif
595	>
596		end subroutine init_FF
335	–
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 345 | 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 383 | 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 397 | 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 405 | 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 449 | 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 466 | 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 476 | 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 493 | 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 502 | 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 526 | 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 567 | 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 589 | 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 610 | 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 625 | 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 645 | 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 681 | 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 718 | 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 792 | 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 805 | Line 1075 | contains
1075		me_i = atid(i)
1076		me_j = atid(j)
1077		#endif
808	–
809	–	if (FF_uses_LJ .and. SIM_uses_LJ) then
810	–
811	–	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
812	–	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
813	–	endif
814	–
815	–	endif
816	–
817	–	if (FF_uses_charges .and. SIM_uses_charges) then
818	–
819	–	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
820	–	call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
821	–	endif
822	–
823	–	endif
824	–
825	–	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
826	–
827	–	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
828	–	call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
829	–	do_pot)
830	–	if (FF_uses_RF .and. SIM_uses_RF) then
831	–	call accumulate_rf(i, j, r, u_l, sw)
832	–	call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
833	–	endif
834	–	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, &
852	<	do_pot)
853	<	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
860	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
861	<	do_pot)
862	<	endif
863	<
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
876	<
877	<	logical, intent(inout) :: do_pot, do_stress
878	<	integer, intent(in) :: i, j
879	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
880	<	real ( kind = dp )                :: r, rc
881	<	real ( kind = dp ), intent(inout) :: d(3), dc(3)
882	<
883	<	logical :: is_EAM_i, is_EAM_j
884	<
885	<	integer :: me_i, me_j
886	<
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
893	<	endif
894	<
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()
978	<
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
980	–
981	–	q_Row = 0.0_dp
982	–	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
1124	–
1125	–	!! This cleans componets of force arrays belonging only to fortran
1380
1381	<	subroutine add_stress_tensor(dpair, fpair)
1128	<
1129	<	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1130	<
1131	<	! because the d vector is the rj - ri vector, and
1132	<	! because fx, fy, fz are the force on atom i, we need a
1133	<	! negative sign here:
1134	<
1135	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1136	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1137	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1138	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1139	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1140	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1141	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1142	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1143	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1144	<
1145	<	virial_Temp = virial_Temp + &
1146	<	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1147	<
1148	<	end subroutine add_stress_tensor
1149	<
1150	<	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(use_RF_c, thisStat)
1155	<	use doForces, ONLY: init_FF
1156	<	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(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)
1165	<
1166	<	use definitions, ONLY: dp
1167	<	use simulation
1168	<	use doForces, ONLY: do_force_loop
1169	<	!! Position array provided by C, dimensioned by getNlocal
1170	<	real ( kind = dp ), dimension(3, nLocal) :: q
1171	<	!! molecular center-of-mass position array
1172	<	real ( kind = dp ), dimension(3, nGroups) :: q_group
1173	<	!! Rotation Matrix for each long range particle in simulation.
1174	<	real( kind = dp), dimension(9, nLocal) :: A
1175	<	!! Unit vectors for dipoles (lab frame)
1176	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1177	<	!! Force array provided by C, dimensioned by getNlocal
1178	<	real ( kind = dp ), dimension(3,nLocal) :: f
1179	<	!! Torsion array provided by C, dimensioned by getNlocal
1180	<	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
1185	<	logical ( kind = 2) :: do_pot_c, do_stress_c
1186	<	integer :: error
1187	<
1188	<	call do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
1189	<	do_pot_c, do_stress_c, error)
1190	<
1191	<	end subroutine doForceloop
1404	>	end subroutine add_stress_tensor
1405	>
1406	>	end module doForces

Diff Legend

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