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

Diff Legend

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