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

Diff Legend

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