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 1650 by gezelter, Tue Oct 26 22:24:52 2004 UTC vs.
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC

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

Diff Legend

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