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 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.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.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 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
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 221 | 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 241 | Line 381 | contains
381		integer :: myStatus
382
383		error = 0
244	–
245	–	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 286 | Line 426 | contains
426		#endif
427		return
428		end subroutine doReadyCheck
289	–
429
430	+
431		subroutine init_FF(use_RF_c, thisStat)
432
433		logical, intent(in) :: use_RF_c
#	Line 302 | 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_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 327 | 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
#	Line 340 | Line 481 | contains
481		FF_uses_Charges = .true.
482		FF_uses_Electrostatics = .true.
483		endif
484	<
484	>
485		call getMatchingElementList(atypes, "is_Dipole", .true., &
486		nMatches, MatchList)
487		if (nMatches .gt. 0) then
#	Line 356 | Line 497 | contains
497		FF_uses_Electrostatics = .true.
498		FF_uses_DirectionalAtoms = .true.
499		endif
500	<
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
#	Line 370 | 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 387 | 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 402 | 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 415 | 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 435 | 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 445 | Line 593 | contains
593		return
594		endif
595		haveNeighborList = .true.
596	<	endif
597	<
596	>	endif
597	>
598		end subroutine init_FF
451	–
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 499 | 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 513 | 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 521 | 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 565 | 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 582 | 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 591 | 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 609 | 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 618 | 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 642 | 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 705 | 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 726 | 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 741 | 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 761 | 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 797 | 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 834 | 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 858 | 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 908 | 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 922 | Line 1082 | contains
1082		me_j = atid(j)
1083		#endif
1084
1085	<	!    write(*,*) i, j, me_i, me_j
1086	<
1087	<	if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1088	<
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	<
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
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
1090
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 (FF_uses_Sticky .and. SIM_uses_sticky) then
1095	>	if (FF_uses_RF .and. SIM_uses_RF) then
1096
1097	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1098	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1099	<	pot, A, f, t, do_pot)
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. &
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	<
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 ( 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 (FF_uses_EAM .and. SIM_uses_EAM) then
1120	<
1121	<	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	<
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
991	<	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
992	<	pot, A, f, t, do_pot)
993	<	endif
994	<
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 999 | 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
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	<
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
1024	<	endif
1025	<
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
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
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)
1243	<	#else
1244	<	!! in the normal loop, the atom numbers are unique
1245	<	unique_id_1 = 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	<	!! We were called with only one atom, so just check the global exclude
1273	<	!! list for this atom
1274	<	if (.not. present(atom2)) then
1275	<	do i = 1, nExcludes_global
1276	<	if (excludesGlobal(i) == unique_id_1) then
1277	<	skip_it = .true.
1278	<	return
1279	<	end if
1280	<	end do
1281	<	return
1282	<	end if
1283	<
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	<	unique_id_2 = AtomColToGlobal(atom2)
1298	>	!! in MPI, we have to look up the unique IDs for each atom
1299	>	unique_id_1 = AtomRowToGlobal(atom1)
1300		#else
1301	<	unique_id_2 = atom2
1301	>	!! in the normal loop, the atom numbers are unique
1302	>	unique_id_1 = atom1
1303		#endif
1304	<
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	<	!! this situation should only arise in MPI simulations
1319	<	if (unique_id_1 == unique_id_2) then
1320	<	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
1318	>	unique_id_2 = AtomColToGlobal(atom2)
1319	>	#else
1320	>	unique_id_2 = atom2
1321		#endif
1322	<
1323	<	!! the rest of these situations can happen in all simulations:
1324	<	do i = 1, nExcludes_global
1325	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1326	<	(excludesGlobal(i) == unique_id_2)) then
1327	<	skip_it = .true.
1328	<	return
1329	<	endif
1330	<	enddo
1331	<
1332	<	do i = 1, nSkipsForAtom(atom1)
1333	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1334	<	skip_it = .true.
1335	<	return
1336	<	endif
1337	<	end do
1338	<
1339	<	return
1340	<	end function skipThisPair
1341	<
1342	<	function FF_UsesDirectionalAtoms() result(doesit)
1343	<	logical :: doesit
1344	<	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1345	<	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1346	<	FF_uses_GayBerne .or. FF_uses_Shapes
1347	<	end function FF_UsesDirectionalAtoms
1348	<
1349	<	function FF_RequiresPrepairCalc() result(doesit)
1350	<	logical :: doesit
1351	<	doesit = FF_uses_EAM
1352	<	end function FF_RequiresPrepairCalc
1353	<
1354	<	function FF_RequiresPostpairCalc() result(doesit)
1355	<	logical :: doesit
1356	<	doesit = FF_uses_RF
1357	<	end function FF_RequiresPostpairCalc
1358	<
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
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_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
1256	–
1257	–	!! 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