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 1706 by gezelter, Thu Nov 4 16:20:28 2004 UTC vs.
Revision 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC

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

Diff Legend

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