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 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2266 by chuckv, Thu Jul 28 22:12:45 2005 UTC

#	Line 45 | Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.10 2005-01-14 20:31:12 gezelter Exp $, $Date: 2005-01-14 20:31:12 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
48	>	!! @version $Id: doForces.F90,v 1.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
49
50
51		module doForces
#	Line 57 | Line 57 | module doForces
57		use neighborLists
58		use lj
59		use sticky
60	<	use dipole_dipole
61	<	use charge_charge
60	>	use electrostatic_module
61		use reaction_field
62		use gb_pair
63		use shapes
#	Line 74 | 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 81 | Line 81 | module doForces
81		logical, save :: haveRlist = .false.
82		logical, save :: haveNeighborList = .false.
83		logical, save :: haveSIMvariables = .false.
84	–	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 101 | 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 112 | 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 125 | 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.
136	<	logical :: is_EAM           = .false.
137	<	logical :: is_Shape         = .false.
138	<	logical :: is_FLARB         = .false.
139	<	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
142	<
147	>
148		contains
149
150	<	subroutine setRlistDF( this_rlist )
151	<
147	<	real(kind=dp) :: this_rlist
148	<
149	<	rlist = this_rlist
150	<	rlistsq = rlist * rlist
151	<
152	<	haveRlist = .true.
153	<
154	<	end subroutine setRlistDF
155	<
156	<	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
171	–
172	–	if (.not. allocated(PropertyMap)) then
173	–	allocate(PropertyMap(nAtypes))
174	–	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)
181	<	PropertyMap(i)%is_LennardJones = thisProperty
182	<
183	<	call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
184	<	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
188	<
189	<	call getElementProperty(atypes, i, "is_Dipole", thisProperty)
190	<	PropertyMap(i)%is_Dipole = thisProperty
204	>	iHash = 0
205	>	myRcut = 0.0_dp
206
207	<	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
208	<	PropertyMap(i)%is_Sticky = 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_GayBerne", thisProperty)
216	<	PropertyMap(i)%is_GayBerne = 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_EAM", thisProperty)
228	<	PropertyMap(i)%is_EAM = 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_Shape", thisProperty)
232	<	PropertyMap(i)%is_Shape = thisProperty
231	>	if (i_is_EAM .and. j_is_EAM) then
232	>	iHash = ior(iHash, EAM_PAIR)
233	>	endif
234
235	<	call getElementProperty(atypes, i, "is_FLARB", thisProperty)
236	<	PropertyMap(i)%is_FLARB = thisProperty
235	>	if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
236	>	if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
237	>	if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
238	>
239	>	if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
240	>	if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
241	>	if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
242	>
243	>
244	>	InteractionMap(i,j)%InteractionHash = iHash
245	>	InteractionMap(j,i)%InteractionHash = iHash
246	>
247	>	end do
248	>
249		end do
250	+	end subroutine createInteractionMap
251
252	<	havePropertyMap = .true.
252	>	! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
253	>	subroutine createRcuts(defaultRList,stat)
254	>	real(kind=dp), intent(in), optional :: defaultRList
255	>	integer :: iMap
256	>	integer :: map_i,map_j
257	>	real(kind=dp) :: thisRCut = 0.0_dp
258	>	real(kind=dp) :: actualCutoff = 0.0_dp
259	>	integer, intent(out) :: stat
260	>	integer :: nAtypes
261	>	integer :: myStatus
262
263	<	end subroutine createPropertyMap
263	>	stat = 0
264	>	if (.not. haveInteractionMap) then
265
266	+	call createInteractionMap(myStatus)
267	+
268	+	if (myStatus .ne. 0) then
269	+	write(default_error, *) 'createInteractionMap failed in doForces!'
270	+	stat = -1
271	+	return
272	+	endif
273	+	endif
274	+
275	+
276	+	nAtypes = getSize(atypes)
277	+	! If we pass a default rcut, set all atypes to that cutoff distance
278	+	if(present(defaultRList)) then
279	+	InteractionMap(:,:)%rList = defaultRList
280	+	InteractionMap(:,:)%rListSq = defaultRList*defaultRList
281	+	haveRlist = .true.
282	+	return
283	+	end if
284	+
285	+	do map_i = 1,nAtypes
286	+	do map_j = map_i,nAtypes
287	+	iMap = InteractionMap(map_i, map_j)%InteractionHash
288	+
289	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
290	+	!            thisRCut = getLJCutOff(map_i,map_j)
291	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
292	+	endif
293	+
294	+	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
295	+	!            thisRCut = getElectrostaticCutOff(map_i,map_j)
296	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
297	+	endif
298	+
299	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300	+	!             thisRCut = getStickyCutOff(map_i,map_j)
301	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
302	+	endif
303	+
304	+	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305	+	!              thisRCut = getStickyPowerCutOff(map_i,map_j)
306	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
307	+	endif
308	+
309	+	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310	+	!              thisRCut = getGayberneCutOff(map_i,map_j)
311	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
312	+	endif
313	+
314	+	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315	+	!              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
317	+	endif
318	+
319	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
320	+	!              thisRCut = getEAMCutOff(map_i,map_j)
321	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
322	+	endif
323	+
324	+	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
325	+	!              thisRCut = getShapeCutOff(map_i,map_j)
326	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
327	+	endif
328	+
329	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
330	+	!              thisRCut = getShapeLJCutOff(map_i,map_j)
331	+	if (thisRcut > actualCutoff) actualCutoff = thisRcut
332	+	endif
333	+	InteractionMap(map_i, map_j)%rList = actualCutoff
334	+	InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335	+	end do
336	+	end do
337	+	haveRlist = .true.
338	+	end subroutine createRcuts
339	+
340	+
341	+	!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
342	+	!!$  subroutine setRlistDF( this_rlist )
343	+	!!$
344	+	!!$   real(kind=dp) :: this_rlist
345	+	!!$
346	+	!!$    rlist = this_rlist
347	+	!!$    rlistsq = rlist * rlist
348	+	!!$
349	+	!!$    haveRlist = .true.
350	+	!!$
351	+	!!$  end subroutine setRlistDF
352	+
353	+
354		subroutine setSimVariables()
355		SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
356		SIM_uses_LennardJones = SimUsesLennardJones()
#	Line 216 | Line 358 | contains
358		SIM_uses_Charges = SimUsesCharges()
359		SIM_uses_Dipoles = SimUsesDipoles()
360		SIM_uses_Sticky = SimUsesSticky()
361	+	SIM_uses_StickyPower = SimUsesStickyPower()
362		SIM_uses_GayBerne = SimUsesGayBerne()
363		SIM_uses_EAM = SimUsesEAM()
364		SIM_uses_Shapes = SimUsesShapes()
#	Line 236 | Line 379 | contains
379		integer :: myStatus
380
381		error = 0
239	–
240	–	if (.not. havePropertyMap) then
382
383	+	if (.not. haveInteractionMap) then
384	+
385		myStatus = 0
386
387	<	call createPropertyMap(myStatus)
387	>	call createInteractionMap(myStatus)
388
389		if (myStatus .ne. 0) then
390	<	write(default_error, *) 'createPropertyMap failed in doForces!'
390	>	write(default_error, *) 'createInteractionMap failed in doForces!'
391		error = -1
392		return
393		endif
#	Line 281 | Line 424 | contains
424		#endif
425		return
426		end subroutine doReadyCheck
284	–
427
428	+
429		subroutine init_FF(use_RF_c, thisStat)
430
431		logical, intent(in) :: use_RF_c
#	Line 297 | Line 440 | contains
440
441		!! Fortran's version of a cast:
442		FF_uses_RF = use_RF_c
443	<
443	>
444		!! init_FF is called *after* all of the atom types have been
445		!! defined in atype_module using the new_atype subroutine.
446		!!
447		!! this will scan through the known atypes and figure out what
448		!! interactions are used by the force field.
449	<
449	>
450		FF_uses_DirectionalAtoms = .false.
451		FF_uses_LennardJones = .false.
452	<	FF_uses_Electrostatic = .false.
452	>	FF_uses_Electrostatics = .false.
453		FF_uses_Charges = .false.
454		FF_uses_Dipoles = .false.
455		FF_uses_Sticky = .false.
456	+	FF_uses_StickyPower = .false.
457		FF_uses_GayBerne = .false.
458		FF_uses_EAM = .false.
459		FF_uses_Shapes = .false.
460		FF_uses_FLARB = .false.
461	<
461	>
462		call getMatchingElementList(atypes, "is_Directional", .true., &
463		nMatches, MatchList)
464		if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
#	Line 322 | Line 466 | contains
466		call getMatchingElementList(atypes, "is_LennardJones", .true., &
467		nMatches, MatchList)
468		if (nMatches .gt. 0) FF_uses_LennardJones = .true.
469	<
469	>
470		call getMatchingElementList(atypes, "is_Electrostatic", .true., &
471		nMatches, MatchList)
472		if (nMatches .gt. 0) then
473	<	FF_uses_Electrostatic = .true.
473	>	FF_uses_Electrostatics = .true.
474		endif
475
476		call getMatchingElementList(atypes, "is_Charge", .true., &
477		nMatches, MatchList)
478		if (nMatches .gt. 0) then
479	<	FF_uses_charges = .true.
480	<	FF_uses_electrostatic = .true.
479	>	FF_uses_Charges = .true.
480	>	FF_uses_Electrostatics = .true.
481		endif
482	<
482	>
483		call getMatchingElementList(atypes, "is_Dipole", .true., &
484		nMatches, MatchList)
485		if (nMatches .gt. 0) then
486	<	FF_uses_dipoles = .true.
487	<	FF_uses_electrostatic = .true.
486	>	FF_uses_Dipoles = .true.
487	>	FF_uses_Electrostatics = .true.
488		FF_uses_DirectionalAtoms = .true.
489		endif
490	<
490	>
491	>	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
492	>	nMatches, MatchList)
493	>	if (nMatches .gt. 0) then
494	>	FF_uses_Quadrupoles = .true.
495	>	FF_uses_Electrostatics = .true.
496	>	FF_uses_DirectionalAtoms = .true.
497	>	endif
498	>
499		call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
500		MatchList)
501		if (nMatches .gt. 0) then
502		FF_uses_Sticky = .true.
503		FF_uses_DirectionalAtoms = .true.
504		endif
505	+
506	+	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
507	+	MatchList)
508	+	if (nMatches .gt. 0) then
509	+	FF_uses_StickyPower = .true.
510	+	FF_uses_DirectionalAtoms = .true.
511	+	endif
512
513		call getMatchingElementList(atypes, "is_GayBerne", .true., &
514		nMatches, MatchList)
#	Line 357 | Line 516 | contains
516		FF_uses_GayBerne = .true.
517		FF_uses_DirectionalAtoms = .true.
518		endif
519	<
519	>
520		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
521		if (nMatches .gt. 0) FF_uses_EAM = .true.
522	<
522	>
523		call getMatchingElementList(atypes, "is_Shape", .true., &
524		nMatches, MatchList)
525		if (nMatches .gt. 0) then
#	Line 374 | Line 533 | contains
533
534		!! Assume sanity (for the sake of argument)
535		haveSaneForceField = .true.
536	<
536	>
537		!! check to make sure the FF_uses_RF setting makes sense
538	<
538	>
539		if (FF_uses_dipoles) then
540		if (FF_uses_RF) then
541		dielect = getDielect()
#	Line 389 | Line 548 | contains
548		haveSaneForceField = .false.
549		return
550		endif
551	<	endif
551	>	endif
552
553		!sticky module does not contain check_sticky_FF anymore
554		!if (FF_uses_sticky) then
#	Line 402 | Line 561 | contains
561		!endif
562
563		if (FF_uses_EAM) then
564	<	call init_EAM_FF(my_status)
564	>	call init_EAM_FF(my_status)
565		if (my_status /= 0) then
566		write(default_error, *) "init_EAM_FF returned a bad status"
567		thisStat = -1
#	Line 422 | Line 581 | contains
581
582		if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583		endif
584	<
584	>
585		if (.not. haveNeighborList) then
586		!! Create neighbor lists
587		call expandNeighborList(nLocal, my_status)
#	Line 432 | Line 591 | contains
591		return
592		endif
593		haveNeighborList = .true.
594	<	endif
595	<
594	>	endif
595	>
596		end subroutine init_FF
438	–
597
598	+
599		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
600		!------------------------------------------------------------->
601		subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
#	Line 486 | Line 645 | contains
645		integer :: localError
646		integer :: propPack_i, propPack_j
647		integer :: loopStart, loopEnd, loop
648	<
648	>	integer :: iMap
649		real(kind=dp) :: listSkin = 1.0
650	<
650	>
651		!! initialize local variables
652	<
652	>
653		#ifdef IS_MPI
654		pot_local = 0.0_dp
655		nAtomsInRow   = getNatomsInRow(plan_atom_row)
#	Line 500 | Line 659 | contains
659		#else
660		natoms = nlocal
661		#endif
662	<
662	>
663		call doReadyCheck(localError)
664		if ( localError .ne. 0 ) then
665		call handleError("do_force_loop", "Not Initialized")
#	Line 508 | Line 667 | contains
667		return
668		end if
669		call zero_work_arrays()
670	<
670	>
671		do_pot = do_pot_c
672		do_stress = do_stress_c
673	<
673	>
674		! Gather all information needed by all force loops:
675	<
675	>
676		#ifdef IS_MPI
677	<
677	>
678		call gather(q, q_Row, plan_atom_row_3d)
679		call gather(q, q_Col, plan_atom_col_3d)
680
681		call gather(q_group, q_group_Row, plan_group_row_3d)
682		call gather(q_group, q_group_Col, plan_group_col_3d)
683	<
683	>
684		if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
685		call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
686		call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
687	<
687	>
688		call gather(A, A_Row, plan_atom_row_rotation)
689		call gather(A, A_Col, plan_atom_col_rotation)
690		endif
691	<
691	>
692		#endif
693	<
693	>
694		!! Begin force loop timing:
695		#ifdef PROFILE
696		call cpu_time(forceTimeInitial)
697		nloops = nloops + 1
698		#endif
699	<
699	>
700		loopEnd = PAIR_LOOP
701		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
702		loopStart = PREPAIR_LOOP
#	Line 552 | Line 711 | contains
711		if (loop .eq. loopStart) then
712		#ifdef IS_MPI
713		call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
714	<	update_nlist)
714	>	update_nlist)
715		#else
716		call checkNeighborList(nGroups, q_group, listSkin, &
717	<	update_nlist)
717	>	update_nlist)
718		#endif
719		endif
720	<
720	>
721		if (update_nlist) then
722		!! save current configuration and construct neighbor list
723		#ifdef IS_MPI
#	Line 569 | Line 728 | contains
728		neighborListSize = size(list)
729		nlist = 0
730		endif
731	<
731	>
732		istart = 1
733		#ifdef IS_MPI
734		iend = nGroupsInRow
#	Line 579 | Line 738 | contains
738		outer: do i = istart, iend
739
740		if (update_nlist) point(i) = nlist + 1
741	<
741	>
742		n_in_i = groupStartRow(i+1) - groupStartRow(i)
743	<
743	>
744		if (update_nlist) then
745		#ifdef IS_MPI
746	+	me_i = atid_row(i)
747		jstart = 1
748		jend = nGroupsInCol
749		#else
750	+	me_i = atid(i)
751		jstart = i+1
752		jend = nGroups
753		#endif
#	Line 596 | Line 757 | contains
757		! make sure group i has neighbors
758		if (jstart .gt. jend) cycle outer
759		endif
760	<
760	>
761		do jnab = jstart, jend
762		if (update_nlist) then
763		j = jnab
#	Line 605 | Line 766 | contains
766		endif
767
768		#ifdef IS_MPI
769	+	me_j = atid_col(j)
770		call get_interatomic_vector(q_group_Row(:,i), &
771		q_group_Col(:,j), d_grp, rgrpsq)
772		#else
773	+	me_j = atid(j)
774		call get_interatomic_vector(q_group(:,i), &
775		q_group(:,j), d_grp, rgrpsq)
776		#endif
777
778	<	if (rgrpsq < rlistsq) then
778	>	if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
779		if (update_nlist) then
780		nlist = nlist + 1
781	<
781	>
782		if (nlist > neighborListSize) then
783		#ifdef IS_MPI
784		call expandNeighborList(nGroupsInRow, listerror)
#	Line 629 | Line 792 | contains
792		end if
793		neighborListSize = size(list)
794		endif
795	<
795	>
796		list(nlist) = j
797		endif
798	<
798	>
799		if (loop .eq. PAIR_LOOP) then
800		vij = 0.0d0
801		fij(1:3) = 0.0d0
802		endif
803	<
803	>
804		call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
805		in_switching_region)
806	<
806	>
807		n_in_j = groupStartCol(j+1) - groupStartCol(j)
808	<
808	>
809		do ia = groupStartRow(i), groupStartRow(i+1)-1
810	<
810	>
811		atom1 = groupListRow(ia)
812	<
812	>
813		inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
814	<
814	>
815		atom2 = groupListCol(jb)
816	<
816	>
817		if (skipThisPair(atom1, atom2)) cycle inner
818
819		if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
#	Line 692 | Line 855 | contains
855		endif
856		enddo inner
857		enddo
858	<
858	>
859		if (loop .eq. PAIR_LOOP) then
860		if (in_switching_region) then
861		swderiv = vij*dswdr/rgrp
862		fij(1) = fij(1) + swderiv*d_grp(1)
863		fij(2) = fij(2) + swderiv*d_grp(2)
864		fij(3) = fij(3) + swderiv*d_grp(3)
865	<
865	>
866		do ia=groupStartRow(i), groupStartRow(i+1)-1
867		atom1=groupListRow(ia)
868		mf = mfactRow(atom1)
#	Line 713 | Line 876 | contains
876		f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
877		#endif
878		enddo
879	<
879	>
880		do jb=groupStartCol(j), groupStartCol(j+1)-1
881		atom2=groupListCol(jb)
882		mf = mfactCol(atom2)
#	Line 728 | Line 891 | contains
891		#endif
892		enddo
893		endif
894	<
894	>
895		if (do_stress) call add_stress_tensor(d_grp, fij)
896		endif
897		end if
898		enddo
899		enddo outer
900	<
900	>
901		if (update_nlist) then
902		#ifdef IS_MPI
903		point(nGroupsInRow + 1) = nlist + 1
#	Line 748 | Line 911 | contains
911		update_nlist = .false.
912		endif
913		endif
914	<
914	>
915		if (loop .eq. PREPAIR_LOOP) then
916		call do_preforce(nlocal, pot)
917		endif
918	<
918	>
919		enddo
920	<
920	>
921		!! Do timing
922		#ifdef PROFILE
923		call cpu_time(forceTimeFinal)
924		forceTime = forceTime + forceTimeFinal - forceTimeInitial
925		#endif
926	<
926	>
927		#ifdef IS_MPI
928		!!distribute forces
929	<
929	>
930		f_temp = 0.0_dp
931		call scatter(f_Row,f_temp,plan_atom_row_3d)
932		do i = 1,nlocal
933		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
934		end do
935	<
935	>
936		f_temp = 0.0_dp
937		call scatter(f_Col,f_temp,plan_atom_col_3d)
938		do i = 1,nlocal
939		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
940		end do
941	<
941	>
942		if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
943		t_temp = 0.0_dp
944		call scatter(t_Row,t_temp,plan_atom_row_3d)
#	Line 784 | Line 947 | contains
947		end do
948		t_temp = 0.0_dp
949		call scatter(t_Col,t_temp,plan_atom_col_3d)
950	<
950	>
951		do i = 1,nlocal
952		t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
953		end do
954		endif
955	<
955	>
956		if (do_pot) then
957		! scatter/gather pot_row into the members of my column
958		call scatter(pot_Row, pot_Temp, plan_atom_row)
959	<
959	>
960		! scatter/gather pot_local into all other procs
961		! add resultant to get total pot
962		do i = 1, nlocal
963		pot_local = pot_local + pot_Temp(i)
964		enddo
965	<
965	>
966		pot_Temp = 0.0_DP
967	<
967	>
968		call scatter(pot_Col, pot_Temp, plan_atom_col)
969		do i = 1, nlocal
970		pot_local = pot_local + pot_Temp(i)
971		enddo
972	<
972	>
973		endif
974		#endif
975	<
975	>
976		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
977	<
977	>
978		if (FF_uses_RF .and. SIM_uses_RF) then
979	<
979	>
980		#ifdef IS_MPI
981		call scatter(rf_Row,rf,plan_atom_row_3d)
982		call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
#	Line 821 | Line 984 | contains
984		rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
985		end do
986		#endif
987	<
987	>
988		do i = 1, nLocal
989	<
989	>
990		rfpot = 0.0_DP
991		#ifdef IS_MPI
992		me_i = atid_row(i)
993		#else
994		me_i = atid(i)
995		#endif
996	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
997
998	<	if (PropertyMap(me_i)%is_Dipole) then
999	<
998	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
999	>
1000		mu_i = getDipoleMoment(me_i)
1001	<
1001	>
1002		!! The reaction field needs to include a self contribution
1003		!! to the field:
1004		call accumulate_self_rf(i, mu_i, eFrame)
#	Line 845 | Line 1009 | contains
1009		pot_local = pot_local + rfpot
1010		#else
1011		pot = pot + rfpot
1012	<
1012	>
1013		#endif
1014	<	endif
1014	>	endif
1015		enddo
1016		endif
1017		endif
1018	<
1019	<
1018	>
1019	>
1020		#ifdef IS_MPI
1021	<
1021	>
1022		if (do_pot) then
1023		pot = pot + pot_local
1024		!! we assume the c code will do the allreduce to get the total potential
1025		!! we could do it right here if we needed to...
1026		endif
1027	<
1027	>
1028		if (do_stress) then
1029		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1030		mpi_comm_world,mpi_err)
1031		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1032		mpi_comm_world,mpi_err)
1033		endif
1034	<
1034	>
1035		#else
1036	<
1036	>
1037		if (do_stress) then
1038		tau = tau_Temp
1039		virial = virial_Temp
1040		endif
1041	<
1041	>
1042		#endif
1043	<
1043	>
1044		end subroutine do_force_loop
1045	<
1045	>
1046		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1047		eFrame, A, f, t, pot, vpair, fpair)
1048
#	Line 895 | Line 1059 | contains
1059		real ( kind = dp ), intent(inout) :: rijsq
1060		real ( kind = dp )                :: r
1061		real ( kind = dp ), intent(inout) :: d(3)
1062	+	real ( kind = dp ) :: ebalance
1063		integer :: me_i, me_j
1064
1065	+	integer :: iMap
1066	+
1067		r = sqrt(rijsq)
1068		vpair = 0.0d0
1069		fpair(1:3) = 0.0d0
#	Line 909 | Line 1076 | contains
1076		me_j = atid(j)
1077		#endif
1078
1079	<	!    write(*,*) i, j, me_i, me_j
913	<
914	<	if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
915	<
916	<	if ( PropertyMap(me_i)%is_LennardJones .and. &
917	<	PropertyMap(me_j)%is_LennardJones ) then
918	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
919	<	endif
920	<
921	<	endif
922	<
923	<	if (FF_uses_charges .and. SIM_uses_charges) then
924	<
925	<	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
926	<	call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
927	<	pot, f, do_pot)
928	<	endif
929	<
930	<	endif
931	<
932	<	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
933	<
934	<	if ( PropertyMap(me_i)%is_Dipole .and. PropertyMap(me_j)%is_Dipole) then
935	<	call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
936	<	pot, eFrame, f, t, do_pot)
937	<	if (FF_uses_RF .and. SIM_uses_RF) then
938	<	call accumulate_rf(i, j, r, eFrame, sw)
939	<	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
940	<	endif
941	<	endif
1079	>	iMap = InteractionMap(me_i, me_j)%InteractionHash
1080
1081	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
1082	+	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1083		endif
1084
1085	<	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1085	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1086	>	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087	>	pot, eFrame, f, t, do_pot)
1088
1089	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1090	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1091	<	pot, A, f, t, do_pot)
1089	>	if (FF_uses_RF .and. SIM_uses_RF) then
1090	>
1091	>	! CHECK ME (RF needs to know about all electrostatic types)
1092	>	call accumulate_rf(i, j, r, eFrame, sw)
1093	>	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094		endif
1095	<
1095	>
1096		endif
1097
1098	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1099	+	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1100	+	pot, A, f, t, do_pot)
1101	+	endif
1102
1103	<	if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1104	<
1105	<	if ( PropertyMap(me_i)%is_GayBerne .and. &
958	<	PropertyMap(me_j)%is_GayBerne) then
959	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
960	<	pot, A, f, t, do_pot)
961	<	endif
962	<
1103	>	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1104	>	call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1105	>	pot, A, f, t, do_pot)
1106		endif
1107	<
1108	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1109	<
1110	<	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
968	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
969	<	do_pot)
970	<	endif
971	<
1107	>
1108	>	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1109	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1110	>	pot, A, f, t, do_pot)
1111		endif
1112	+
1113	+	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1114	+	!      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1115	+	!           pot, A, f, t, do_pot)
1116	+	endif
1117
1118	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1119	+	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1120	+	do_pot)
1121	+	endif
1122
1123	<	!    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1123	>	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
1124	>	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1125	>	pot, A, f, t, do_pot)
1126	>	endif
1127
1128	<	if (FF_uses_Shapes .and. SIM_uses_Shapes) then
1129	<	if ( PropertyMap(me_i)%is_Shape .and. &
1130	<	PropertyMap(me_j)%is_Shape ) then
980	<	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
981	<	pot, A, f, t, do_pot)
982	<	endif
983	<
1128	>	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
1129	>	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1130	>	pot, A, f, t, do_pot)
1131		endif
1132
1133		end subroutine do_pair
#	Line 988 | Line 1135 | contains
1135		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1136		do_pot, do_stress, eFrame, A, f, t, pot)
1137
1138	<	real( kind = dp ) :: pot, sw
1139	<	real( kind = dp ), dimension(9,nLocal) :: eFrame
1140	<	real (kind=dp), dimension(9,nLocal) :: A
1141	<	real (kind=dp), dimension(3,nLocal) :: f
1142	<	real (kind=dp), dimension(3,nLocal) :: t
996	<
997	<	logical, intent(inout) :: do_pot, do_stress
998	<	integer, intent(in) :: i, j
999	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1000	<	real ( kind = dp )                :: r, rc
1001	<	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1002	<
1003	<	logical :: is_EAM_i, is_EAM_j
1004	<
1005	<	integer :: me_i, me_j
1006	<
1138	>	real( kind = dp ) :: pot, sw
1139	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1140	>	real (kind=dp), dimension(9,nLocal) :: A
1141	>	real (kind=dp), dimension(3,nLocal) :: f
1142	>	real (kind=dp), dimension(3,nLocal) :: t
1143
1144	<	r = sqrt(rijsq)
1145	<	if (SIM_uses_molecular_cutoffs) then
1146	<	rc = sqrt(rcijsq)
1147	<	else
1148	<	rc = r
1013	<	endif
1014	<
1144	>	logical, intent(inout) :: do_pot, do_stress
1145	>	integer, intent(in) :: i, j
1146	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1147	>	real ( kind = dp )                :: r, rc
1148	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1149
1150	+	integer :: me_i, me_j, iMap
1151	+
1152		#ifdef IS_MPI
1153	<	me_i = atid_row(i)
1154	<	me_j = atid_col(j)
1153	>	me_i = atid_row(i)
1154	>	me_j = atid_col(j)
1155		#else
1156	<	me_i = atid(i)
1157	<	me_j = atid(j)
1156	>	me_i = atid(i)
1157	>	me_j = atid(j)
1158		#endif
1023	–
1024	–	if (FF_uses_EAM .and. SIM_uses_EAM) then
1025	–
1026	–	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1027	–	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1028	–
1029	–	endif
1030	–
1031	–	end subroutine do_prepair
1032	–
1033	–
1034	–	subroutine do_preforce(nlocal,pot)
1035	–	integer :: nlocal
1036	–	real( kind = dp ) :: pot
1037	–
1038	–	if (FF_uses_EAM .and. SIM_uses_EAM) then
1039	–	call calc_EAM_preforce_Frho(nlocal,pot)
1040	–	endif
1041	–
1042	–
1043	–	end subroutine do_preforce
1044	–
1045	–
1046	–	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1047	–
1048	–	real (kind = dp), dimension(3) :: q_i
1049	–	real (kind = dp), dimension(3) :: q_j
1050	–	real ( kind = dp ), intent(out) :: r_sq
1051	–	real( kind = dp ) :: d(3), scaled(3)
1052	–	integer i
1053	–
1054	–	d(1:3) = q_j(1:3) - q_i(1:3)
1055	–
1056	–	! Wrap back into periodic box if necessary
1057	–	if ( SIM_uses_PBC ) then
1058	–
1059	–	if( .not.boxIsOrthorhombic ) then
1060	–	! calc the scaled coordinates.
1061	–
1062	–	scaled = matmul(HmatInv, d)
1063	–
1064	–	! wrap the scaled coordinates
1065	–
1066	–	scaled = scaled  - anint(scaled)
1067	–
1068	–
1069	–	! calc the wrapped real coordinates from the wrapped scaled
1070	–	! coordinates
1071	–
1072	–	d = matmul(Hmat,scaled)
1073	–
1074	–	else
1075	–	! calc the scaled coordinates.
1076	–
1077	–	do i = 1, 3
1078	–	scaled(i) = d(i) * HmatInv(i,i)
1079	–
1080	–	! wrap the scaled coordinates
1081	–
1082	–	scaled(i) = scaled(i) - anint(scaled(i))
1083	–
1084	–	! calc the wrapped real coordinates from the wrapped scaled
1085	–	! coordinates
1086	–
1087	–	d(i) = scaled(i)*Hmat(i,i)
1088	–	enddo
1089	–	endif
1090	–
1091	–	endif
1092	–
1093	–	r_sq = dot_product(d,d)
1094	–
1095	–	end subroutine get_interatomic_vector
1096	–
1097	–	subroutine zero_work_arrays()
1098	–
1099	–	#ifdef IS_MPI
1100	–
1101	–	q_Row = 0.0_dp
1102	–	q_Col = 0.0_dp
1159
1160	<	q_group_Row = 0.0_dp
1161	<	q_group_Col = 0.0_dp
1162	<
1163	<	eFrame_Row = 0.0_dp
1164	<	eFrame_Col = 0.0_dp
1165	<
1166	<	A_Row = 0.0_dp
1167	<	A_Col = 0.0_dp
1168	<
1169	<	f_Row = 0.0_dp
1170	<	f_Col = 0.0_dp
1171	<	f_Temp = 0.0_dp
1172	<
1173	<	t_Row = 0.0_dp
1174	<	t_Col = 0.0_dp
1175	<	t_Temp = 0.0_dp
1176	<
1177	<	pot_Row = 0.0_dp
1178	<	pot_Col = 0.0_dp
1179	<	pot_Temp = 0.0_dp
1180	<
1181	<	rf_Row = 0.0_dp
1182	<	rf_Col = 0.0_dp
1183	<	rf_Temp = 0.0_dp
1184	<
1185	<	#endif
1186	<
1187	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1188	<	call clean_EAM()
1189	<	endif
1190	<
1191	<	rf = 0.0_dp
1192	<	tau_Temp = 0.0_dp
1193	<	virial_Temp = 0.0_dp
1194	<	end subroutine zero_work_arrays
1195	<
1196	<	function skipThisPair(atom1, atom2) result(skip_it)
1197	<	integer, intent(in) :: atom1
1198	<	integer, intent(in), optional :: atom2
1199	<	logical :: skip_it
1200	<	integer :: unique_id_1, unique_id_2
1201	<	integer :: me_i,me_j
1202	<	integer :: i
1203	<
1204	<	skip_it = .false.
1205	<
1206	<	!! there are a number of reasons to skip a pair or a particle
1207	<	!! mostly we do this to exclude atoms who are involved in short
1208	<	!! range interactions (bonds, bends, torsions), but we also need
1209	<	!! to exclude some overcounted interactions that result from
1210	<	!! the parallel decomposition
1211	<
1160	>	iMap = InteractionMap(me_i, me_j)%InteractionHash
1161	>
1162	>	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1163	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1164	>	endif
1165	>
1166	>	end subroutine do_prepair
1167	>
1168	>
1169	>	subroutine do_preforce(nlocal,pot)
1170	>	integer :: nlocal
1171	>	real( kind = dp ) :: pot
1172	>
1173	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1174	>	call calc_EAM_preforce_Frho(nlocal,pot)
1175	>	endif
1176	>
1177	>
1178	>	end subroutine do_preforce
1179	>
1180	>
1181	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1182	>
1183	>	real (kind = dp), dimension(3) :: q_i
1184	>	real (kind = dp), dimension(3) :: q_j
1185	>	real ( kind = dp ), intent(out) :: r_sq
1186	>	real( kind = dp ) :: d(3), scaled(3)
1187	>	integer i
1188	>
1189	>	d(1:3) = q_j(1:3) - q_i(1:3)
1190	>
1191	>	! Wrap back into periodic box if necessary
1192	>	if ( SIM_uses_PBC ) then
1193	>
1194	>	if( .not.boxIsOrthorhombic ) then
1195	>	! calc the scaled coordinates.
1196	>
1197	>	scaled = matmul(HmatInv, d)
1198	>
1199	>	! wrap the scaled coordinates
1200	>
1201	>	scaled = scaled  - anint(scaled)
1202	>
1203	>
1204	>	! calc the wrapped real coordinates from the wrapped scaled
1205	>	! coordinates
1206	>
1207	>	d = matmul(Hmat,scaled)
1208	>
1209	>	else
1210	>	! calc the scaled coordinates.
1211	>
1212	>	do i = 1, 3
1213	>	scaled(i) = d(i) * HmatInv(i,i)
1214	>
1215	>	! wrap the scaled coordinates
1216	>
1217	>	scaled(i) = scaled(i) - anint(scaled(i))
1218	>
1219	>	! calc the wrapped real coordinates from the wrapped scaled
1220	>	! coordinates
1221	>
1222	>	d(i) = scaled(i)*Hmat(i,i)
1223	>	enddo
1224	>	endif
1225	>
1226	>	endif
1227	>
1228	>	r_sq = dot_product(d,d)
1229	>
1230	>	end subroutine get_interatomic_vector
1231	>
1232	>	subroutine zero_work_arrays()
1233	>
1234		#ifdef IS_MPI
1235	<	!! in MPI, we have to look up the unique IDs for each atom
1236	<	unique_id_1 = AtomRowToGlobal(atom1)
1235	>
1236	>	q_Row = 0.0_dp
1237	>	q_Col = 0.0_dp
1238	>
1239	>	q_group_Row = 0.0_dp
1240	>	q_group_Col = 0.0_dp
1241	>
1242	>	eFrame_Row = 0.0_dp
1243	>	eFrame_Col = 0.0_dp
1244	>
1245	>	A_Row = 0.0_dp
1246	>	A_Col = 0.0_dp
1247	>
1248	>	f_Row = 0.0_dp
1249	>	f_Col = 0.0_dp
1250	>	f_Temp = 0.0_dp
1251	>
1252	>	t_Row = 0.0_dp
1253	>	t_Col = 0.0_dp
1254	>	t_Temp = 0.0_dp
1255	>
1256	>	pot_Row = 0.0_dp
1257	>	pot_Col = 0.0_dp
1258	>	pot_Temp = 0.0_dp
1259	>
1260	>	rf_Row = 0.0_dp
1261	>	rf_Col = 0.0_dp
1262	>	rf_Temp = 0.0_dp
1263	>
1264	>	#endif
1265	>
1266	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1267	>	call clean_EAM()
1268	>	endif
1269	>
1270	>	rf = 0.0_dp
1271	>	tau_Temp = 0.0_dp
1272	>	virial_Temp = 0.0_dp
1273	>	end subroutine zero_work_arrays
1274	>
1275	>	function skipThisPair(atom1, atom2) result(skip_it)
1276	>	integer, intent(in) :: atom1
1277	>	integer, intent(in), optional :: atom2
1278	>	logical :: skip_it
1279	>	integer :: unique_id_1, unique_id_2
1280	>	integer :: me_i,me_j
1281	>	integer :: i
1282	>
1283	>	skip_it = .false.
1284	>
1285	>	!! there are a number of reasons to skip a pair or a particle
1286	>	!! mostly we do this to exclude atoms who are involved in short
1287	>	!! range interactions (bonds, bends, torsions), but we also need
1288	>	!! to exclude some overcounted interactions that result from
1289	>	!! the parallel decomposition
1290	>
1291	>	#ifdef IS_MPI
1292	>	!! in MPI, we have to look up the unique IDs for each atom
1293	>	unique_id_1 = AtomRowToGlobal(atom1)
1294		#else
1295	<	!! in the normal loop, the atom numbers are unique
1296	<	unique_id_1 = atom1
1295	>	!! in the normal loop, the atom numbers are unique
1296	>	unique_id_1 = atom1
1297		#endif
1298	<
1299	<	!! We were called with only one atom, so just check the global exclude
1300	<	!! list for this atom
1301	<	if (.not. present(atom2)) then
1302	<	do i = 1, nExcludes_global
1303	<	if (excludesGlobal(i) == unique_id_1) then
1304	<	skip_it = .true.
1305	<	return
1306	<	end if
1307	<	end do
1308	<	return
1309	<	end if
1310	<
1298	>
1299	>	!! We were called with only one atom, so just check the global exclude
1300	>	!! list for this atom
1301	>	if (.not. present(atom2)) then
1302	>	do i = 1, nExcludes_global
1303	>	if (excludesGlobal(i) == unique_id_1) then
1304	>	skip_it = .true.
1305	>	return
1306	>	end if
1307	>	end do
1308	>	return
1309	>	end if
1310	>
1311		#ifdef IS_MPI
1312	<	unique_id_2 = AtomColToGlobal(atom2)
1312	>	unique_id_2 = AtomColToGlobal(atom2)
1313		#else
1314	<	unique_id_2 = atom2
1314	>	unique_id_2 = atom2
1315		#endif
1316	<
1316	>
1317		#ifdef IS_MPI
1318	<	!! this situation should only arise in MPI simulations
1319	<	if (unique_id_1 == unique_id_2) then
1320	<	skip_it = .true.
1321	<	return
1322	<	end if
1323	<
1324	<	!! this prevents us from doing the pair on multiple processors
1325	<	if (unique_id_1 < unique_id_2) then
1326	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1327	<	skip_it = .true.
1328	<	return
1329	<	endif
1330	<	else
1331	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1332	<	skip_it = .true.
1333	<	return
1334	<	endif
1335	<	endif
1318	>	!! this situation should only arise in MPI simulations
1319	>	if (unique_id_1 == unique_id_2) then
1320	>	skip_it = .true.
1321	>	return
1322	>	end if
1323	>
1324	>	!! this prevents us from doing the pair on multiple processors
1325	>	if (unique_id_1 < unique_id_2) then
1326	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1327	>	skip_it = .true.
1328	>	return
1329	>	endif
1330	>	else
1331	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1332	>	skip_it = .true.
1333	>	return
1334	>	endif
1335	>	endif
1336		#endif
1337	<
1338	<	!! the rest of these situations can happen in all simulations:
1339	<	do i = 1, nExcludes_global
1340	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1341	<	(excludesGlobal(i) == unique_id_2)) then
1342	<	skip_it = .true.
1343	<	return
1344	<	endif
1345	<	enddo
1346	<
1347	<	do i = 1, nSkipsForAtom(atom1)
1348	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1349	<	skip_it = .true.
1350	<	return
1351	<	endif
1352	<	end do
1353	<
1354	<	return
1355	<	end function skipThisPair
1356	<
1357	<	function FF_UsesDirectionalAtoms() result(doesit)
1358	<	logical :: doesit
1359	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360	<	FF_uses_Sticky .or. FF_uses_GayBerne .or. FF_uses_Shapes
1361	<	end function FF_UsesDirectionalAtoms
1362	<
1363	<	function FF_RequiresPrepairCalc() result(doesit)
1364	<	logical :: doesit
1365	<	doesit = FF_uses_EAM
1366	<	end function FF_RequiresPrepairCalc
1367	<
1368	<	function FF_RequiresPostpairCalc() result(doesit)
1369	<	logical :: doesit
1370	<	doesit = FF_uses_RF
1371	<	end function FF_RequiresPostpairCalc
1372	<
1337	>
1338	>	!! the rest of these situations can happen in all simulations:
1339	>	do i = 1, nExcludes_global
1340	>	if ((excludesGlobal(i) == unique_id_1) .or. &
1341	>	(excludesGlobal(i) == unique_id_2)) then
1342	>	skip_it = .true.
1343	>	return
1344	>	endif
1345	>	enddo
1346	>
1347	>	do i = 1, nSkipsForAtom(atom1)
1348	>	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1349	>	skip_it = .true.
1350	>	return
1351	>	endif
1352	>	end do
1353	>
1354	>	return
1355	>	end function skipThisPair
1356	>
1357	>	function FF_UsesDirectionalAtoms() result(doesit)
1358	>	logical :: doesit
1359	>	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360	>	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1361	>	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1362	>	end function FF_UsesDirectionalAtoms
1363	>
1364	>	function FF_RequiresPrepairCalc() result(doesit)
1365	>	logical :: doesit
1366	>	doesit = FF_uses_EAM
1367	>	end function FF_RequiresPrepairCalc
1368	>
1369	>	function FF_RequiresPostpairCalc() result(doesit)
1370	>	logical :: doesit
1371	>	doesit = FF_uses_RF
1372	>	end function FF_RequiresPostpairCalc
1373	>
1374		#ifdef PROFILE
1375	<	function getforcetime() result(totalforcetime)
1376	<	real(kind=dp) :: totalforcetime
1377	<	totalforcetime = forcetime
1378	<	end function getforcetime
1375	>	function getforcetime() result(totalforcetime)
1376	>	real(kind=dp) :: totalforcetime
1377	>	totalforcetime = forcetime
1378	>	end function getforcetime
1379		#endif
1244	–
1245	–	!! This cleans componets of force arrays belonging only to fortran
1380
1381	<	subroutine add_stress_tensor(dpair, fpair)
1382	<
1383	<	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1384	<
1385	<	! because the d vector is the rj - ri vector, and
1386	<	! because fx, fy, fz are the force on atom i, we need a
1387	<	! negative sign here:
1388	<
1389	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1390	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1391	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1392	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1393	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1394	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1395	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1396	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1397	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1398	<
1399	<	virial_Temp = virial_Temp + &
1400	<	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1401	<
1402	<	end subroutine add_stress_tensor
1403	<
1381	>	!! This cleans componets of force arrays belonging only to fortran
1382	>
1383	>	subroutine add_stress_tensor(dpair, fpair)
1384	>
1385	>	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1386	>
1387	>	! because the d vector is the rj - ri vector, and
1388	>	! because fx, fy, fz are the force on atom i, we need a
1389	>	! negative sign here:
1390	>
1391	>	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1392	>	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1393	>	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1394	>	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1395	>	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1396	>	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1397	>	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1398	>	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1399	>	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1400	>
1401	>	virial_Temp = virial_Temp + &
1402	>	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1403	>
1404	>	end subroutine add_stress_tensor
1405	>
1406		end module doForces

Diff Legend

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