ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-3.0/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2085 by gezelter, Tue Mar 8 21:05:46 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.11 2005-03-08 21:05:46 gezelter Exp $, $Date: 2005-03-08 21:05:46 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
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 73 | 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 80 | Line 81 | module doForces
81		logical, save :: haveRlist = .false.
82		logical, save :: haveNeighborList = .false.
83		logical, save :: haveSIMvariables = .false.
83	–	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_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
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 103 | Line 105 | module doForces
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 113 | 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 126 | 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_Quadrupole    = .false.
145	<	logical :: is_Sticky        = .false.
137	<	logical :: is_GayBerne      = .false.
138	<	logical :: is_EAM           = .false.
139	<	logical :: is_Shape         = .false.
140	<	logical :: is_FLARB         = .false.
141	<	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
144	<
147	>
148		contains
149
150	<	subroutine setRlistDF( this_rlist )
151	<
149	<	real(kind=dp) :: this_rlist
150	<
151	<	rlist = this_rlist
152	<	rlistsq = rlist * rlist
153	<
154	<	haveRlist = .true.
155	<
156	<	end subroutine setRlistDF
157	<
158	<	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
173	–
174	–	if (.not. allocated(PropertyMap)) then
175	–	allocate(PropertyMap(nAtypes))
176	–	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)
183	<	PropertyMap(i)%is_LennardJones = thisProperty
184	<
185	<	call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
186	<	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
190	<
191	<	call getElementProperty(atypes, i, "is_Dipole", thisProperty)
192	<	PropertyMap(i)%is_Dipole = thisProperty
204	>	iHash = 0
205	>	myRcut = 0.0_dp
206
207	<	call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
208	<	PropertyMap(i)%is_Quadrupole = 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_Sticky", thisProperty)
216	<	PropertyMap(i)%is_Sticky = 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_GayBerne", thisProperty)
228	<	PropertyMap(i)%is_GayBerne = 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_EAM", thisProperty)
232	<	PropertyMap(i)%is_EAM = 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_Shape", thisProperty)
236	<	PropertyMap(i)%is_Shape = 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	<	call getElementProperty(atypes, i, "is_FLARB", thisProperty)
240	<	PropertyMap(i)%is_FLARB = thisProperty
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 221 | 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 241 | Line 379 | contains
379		integer :: myStatus
380
381		error = 0
244	–
245	–	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 286 | Line 424 | contains
424		#endif
425		return
426		end subroutine doReadyCheck
289	–
427
428	+
429		subroutine init_FF(use_RF_c, thisStat)
430
431		logical, intent(in) :: use_RF_c
#	Line 302 | 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_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 327 | 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
#	Line 340 | Line 479 | contains
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
#	Line 356 | Line 495 | contains
495		FF_uses_Electrostatics = .true.
496		FF_uses_DirectionalAtoms = .true.
497		endif
498	<
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 370 | 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 387 | 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 402 | 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 415 | 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 435 | 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 445 | Line 591 | contains
591		return
592		endif
593		haveNeighborList = .true.
594	<	endif
595	<
594	>	endif
595	>
596		end subroutine init_FF
451	–
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 499 | 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 513 | 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 521 | 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 565 | 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 582 | 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 592 | 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 609 | 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 618 | 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 642 | 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 705 | 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 726 | 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 741 | 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 761 | 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 797 | 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 834 | 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 858 | 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 908 | 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 922 | Line 1076 | contains
1076		me_j = atid(j)
1077		#endif
1078
1079	<	!    write(*,*) i, j, me_i, me_j
1080	<
1081	<	if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1082	<
929	<	if ( PropertyMap(me_i)%is_LennardJones .and. &
930	<	PropertyMap(me_j)%is_LennardJones ) then
931	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
932	<	endif
933	<
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
935	–
936	–	if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
937	–
938	–	if (PropertyMap(me_i)%is_Electrostatic .and. &
939	–	PropertyMap(me_j)%is_Electrostatic) then
940	–	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
941	–	pot, eFrame, f, t, do_pot)
942	–	endif
943	–
944	–	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
945	–	if ( PropertyMap(me_i)%is_Dipole .and. &
946	–	PropertyMap(me_j)%is_Dipole) then
947	–	if (FF_uses_RF .and. SIM_uses_RF) then
948	–	call accumulate_rf(i, j, r, eFrame, sw)
949	–	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
950	–	endif
951	–	endif
952	–	endif
953	–	endif
1084
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 (FF_uses_Sticky .and. SIM_uses_sticky) then
1089	>	if (FF_uses_RF .and. SIM_uses_RF) then
1090
1091	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1092	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093	<	pot, A, f, t, do_pot)
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. &
969	<	PropertyMap(me_j)%is_GayBerne) then
970	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
971	<	pot, A, f, t, do_pot)
972	<	endif
973	<
1103	>	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1104	>	call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1105	>	pot, A, f, t, do_pot)
1106		endif
1107	+
1108	+	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1109	+	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1110	+	pot, A, f, t, do_pot)
1111	+	endif
1112
1113	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1114	<
1115	<	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
979	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
980	<	do_pot)
981	<	endif
982	<
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
991	<	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
992	<	pot, A, f, t, do_pot)
993	<	endif
994	<
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 999 | 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
1007	<
1008	<	logical, intent(inout) :: do_pot, do_stress
1009	<	integer, intent(in) :: i, j
1010	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1011	<	real ( kind = dp )                :: r, rc
1012	<	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1013	<
1014	<	logical :: is_EAM_i, is_EAM_j
1015	<
1016	<	integer :: me_i, me_j
1017	<
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
1024	<	endif
1025	<
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
1034	–
1035	–	if (FF_uses_EAM .and. SIM_uses_EAM) then
1036	–
1037	–	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1038	–	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1039	–
1040	–	endif
1041	–
1042	–	end subroutine do_prepair
1043	–
1044	–
1045	–	subroutine do_preforce(nlocal,pot)
1046	–	integer :: nlocal
1047	–	real( kind = dp ) :: pot
1048	–
1049	–	if (FF_uses_EAM .and. SIM_uses_EAM) then
1050	–	call calc_EAM_preforce_Frho(nlocal,pot)
1051	–	endif
1052	–
1053	–
1054	–	end subroutine do_preforce
1055	–
1056	–
1057	–	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1058	–
1059	–	real (kind = dp), dimension(3) :: q_i
1060	–	real (kind = dp), dimension(3) :: q_j
1061	–	real ( kind = dp ), intent(out) :: r_sq
1062	–	real( kind = dp ) :: d(3), scaled(3)
1063	–	integer i
1064	–
1065	–	d(1:3) = q_j(1:3) - q_i(1:3)
1066	–
1067	–	! Wrap back into periodic box if necessary
1068	–	if ( SIM_uses_PBC ) then
1069	–
1070	–	if( .not.boxIsOrthorhombic ) then
1071	–	! calc the scaled coordinates.
1072	–
1073	–	scaled = matmul(HmatInv, d)
1074	–
1075	–	! wrap the scaled coordinates
1076	–
1077	–	scaled = scaled  - anint(scaled)
1078	–
1079	–
1080	–	! calc the wrapped real coordinates from the wrapped scaled
1081	–	! coordinates
1082	–
1083	–	d = matmul(Hmat,scaled)
1084	–
1085	–	else
1086	–	! calc the scaled coordinates.
1087	–
1088	–	do i = 1, 3
1089	–	scaled(i) = d(i) * HmatInv(i,i)
1090	–
1091	–	! wrap the scaled coordinates
1092	–
1093	–	scaled(i) = scaled(i) - anint(scaled(i))
1094	–
1095	–	! calc the wrapped real coordinates from the wrapped scaled
1096	–	! coordinates
1097	–
1098	–	d(i) = scaled(i)*Hmat(i,i)
1099	–	enddo
1100	–	endif
1101	–
1102	–	endif
1103	–
1104	–	r_sq = dot_product(d,d)
1105	–
1106	–	end subroutine get_interatomic_vector
1107	–
1108	–	subroutine zero_work_arrays()
1109	–
1110	–	#ifdef IS_MPI
1111	–
1112	–	q_Row = 0.0_dp
1113	–	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)
1237	<	#else
1238	<	!! in the normal loop, the atom numbers are unique
1239	<	unique_id_1 = 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	<	!! We were called with only one atom, so just check the global exclude
1267	<	!! list for this atom
1268	<	if (.not. present(atom2)) then
1269	<	do i = 1, nExcludes_global
1270	<	if (excludesGlobal(i) == unique_id_1) then
1271	<	skip_it = .true.
1272	<	return
1273	<	end if
1274	<	end do
1275	<	return
1276	<	end if
1277	<
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	<	unique_id_2 = AtomColToGlobal(atom2)
1292	>	!! in MPI, we have to look up the unique IDs for each atom
1293	>	unique_id_1 = AtomRowToGlobal(atom1)
1294		#else
1295	<	unique_id_2 = atom2
1295	>	!! in the normal loop, the atom numbers are unique
1296	>	unique_id_1 = atom1
1297		#endif
1298	<
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	<	!! this situation should only arise in MPI simulations
1313	<	if (unique_id_1 == unique_id_2) then
1314	<	skip_it = .true.
1197	<	return
1198	<	end if
1199	<
1200	<	!! this prevents us from doing the pair on multiple processors
1201	<	if (unique_id_1 < unique_id_2) then
1202	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1203	<	skip_it = .true.
1204	<	return
1205	<	endif
1206	<	else
1207	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1208	<	skip_it = .true.
1209	<	return
1210	<	endif
1211	<	endif
1312	>	unique_id_2 = AtomColToGlobal(atom2)
1313	>	#else
1314	>	unique_id_2 = atom2
1315		#endif
1316	<
1317	<	!! the rest of these situations can happen in all simulations:
1318	<	do i = 1, nExcludes_global
1319	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1320	<	(excludesGlobal(i) == unique_id_2)) then
1321	<	skip_it = .true.
1322	<	return
1323	<	endif
1324	<	enddo
1325	<
1326	<	do i = 1, nSkipsForAtom(atom1)
1327	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1328	<	skip_it = .true.
1329	<	return
1330	<	endif
1331	<	end do
1332	<
1333	<	return
1334	<	end function skipThisPair
1335	<
1336	<	function FF_UsesDirectionalAtoms() result(doesit)
1337	<	logical :: doesit
1338	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1339	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1340	<	FF_uses_GayBerne .or. FF_uses_Shapes
1341	<	end function FF_UsesDirectionalAtoms
1342	<
1343	<	function FF_RequiresPrepairCalc() result(doesit)
1344	<	logical :: doesit
1345	<	doesit = FF_uses_EAM
1346	<	end function FF_RequiresPrepairCalc
1347	<
1348	<	function FF_RequiresPostpairCalc() result(doesit)
1349	<	logical :: doesit
1350	<	doesit = FF_uses_RF
1351	<	end function FF_RequiresPostpairCalc
1352	<
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
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_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
1256	–
1257	–	!! 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