ViewVC Help

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

root/group/trunk/OOPSE/libmdtools/do_Forces.F90

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 482 by chuckv, Tue Apr 8 22:38:43 2003 UTC vs.
Revision 1144 by tim, Sat May 1 18:52:38 2004 UTC

#	Line 4 | Line 4
4
5		!! @author Charles F. Vardeman II
6		!! @author Matthew Meineke
7	<	!! @version $Id: do_Forces.F90,v 1.13 2003-04-08 22:38:43 chuckv Exp $, $Date: 2003-04-08 22:38:43 $, $Name: not supported by cvs2svn $, $Revision: 1.13 $
7	>	!! @version $Id: do_Forces.F90,v 1.53 2004-05-01 18:52:38 tim Exp $, $Date: 2004-05-01 18:52:38 $, $Name: not supported by cvs2svn $, $Revision: 1.53 $
8
9		module do_Forces
10		use force_globals
#	Line 15 | Line 15 | module do_Forces
15		use lj
16		use sticky_pair
17		use dipole_dipole
18	+	use charge_charge
19		use reaction_field
20		use gb_pair
21	+	use vector_class
22	+	use eam
23	+	use status
24		#ifdef IS_MPI
25		use mpiSimulation
26		#endif
#	Line 27 | Line 31 | module do_Forces
31		#define __FORTRAN90
32		#include "fForceField.h"
33
34	<	logical, save :: do_forces_initialized = .false.
34	>	logical, save :: haveRlist = .false.
35	>	logical, save :: haveNeighborList = .false.
36	>	logical, save :: havePolicies = .false.
37	>	logical, save :: haveSIMvariables = .false.
38	>	logical, save :: havePropertyMap = .false.
39	>	logical, save :: haveSaneForceField = .false.
40		logical, save :: FF_uses_LJ
41		logical, save :: FF_uses_sticky
42	+	logical, save :: FF_uses_charges
43		logical, save :: FF_uses_dipoles
44		logical, save :: FF_uses_RF
45		logical, save :: FF_uses_GB
46		logical, save :: FF_uses_EAM
47	+	logical, save :: SIM_uses_LJ
48	+	logical, save :: SIM_uses_sticky
49	+	logical, save :: SIM_uses_charges
50	+	logical, save :: SIM_uses_dipoles
51	+	logical, save :: SIM_uses_RF
52	+	logical, save :: SIM_uses_GB
53	+	logical, save :: SIM_uses_EAM
54	+	logical, save :: SIM_requires_postpair_calc
55	+	logical, save :: SIM_requires_prepair_calc
56	+	logical, save :: SIM_uses_directional_atoms
57	+	logical, save :: SIM_uses_PBC
58	+	logical, save :: SIM_uses_molecular_cutoffs
59
60	+	real(kind=dp), save :: rlist, rlistsq
61	+
62		public :: init_FF
63		public :: do_force_loop
64	+	public :: setRlistDF
65
66	+
67	+	#ifdef PROFILE
68	+	public :: getforcetime
69	+	real, save :: forceTime = 0
70	+	real :: forceTimeInitial, forceTimeFinal
71	+	integer :: nLoops
72	+	#endif
73	+
74	+	type :: Properties
75	+	logical :: is_lj     = .false.
76	+	logical :: is_sticky = .false.
77	+	logical :: is_dp     = .false.
78	+	logical :: is_gb     = .false.
79	+	logical :: is_eam    = .false.
80	+	logical :: is_charge = .false.
81	+	real(kind=DP) :: charge = 0.0_DP
82	+	real(kind=DP) :: dipole_moment = 0.0_DP
83	+	end type Properties
84	+
85	+	type(Properties), dimension(:),allocatable :: PropertyMap
86	+
87		contains
88	+
89	+	subroutine setRlistDF( this_rlist )
90	+
91	+	real(kind=dp) :: this_rlist
92	+
93	+	rlist = this_rlist
94	+	rlistsq = rlist * rlist
95	+
96	+	haveRlist = .true.
97	+
98	+	end subroutine setRlistDF
99	+
100	+	subroutine createPropertyMap(status)
101	+	integer :: nAtypes
102	+	integer :: status
103	+	integer :: i
104	+	logical :: thisProperty
105	+	real (kind=DP) :: thisDPproperty
106	+
107	+	status = 0
108	+
109	+	nAtypes = getSize(atypes)
110	+
111	+	if (nAtypes == 0) then
112	+	status = -1
113	+	return
114	+	end if
115	+
116	+	if (.not. allocated(PropertyMap)) then
117	+	allocate(PropertyMap(nAtypes))
118	+	endif
119	+
120	+	do i = 1, nAtypes
121	+	call getElementProperty(atypes, i, "is_LJ", thisProperty)
122	+	PropertyMap(i)%is_LJ = thisProperty
123	+
124	+	call getElementProperty(atypes, i, "is_Charge", thisProperty)
125	+	PropertyMap(i)%is_Charge = thisProperty
126	+
127	+	if (thisProperty) then
128	+	call getElementProperty(atypes, i, "charge", thisDPproperty)
129	+	PropertyMap(i)%charge = thisDPproperty
130	+	endif
131	+
132	+	call getElementProperty(atypes, i, "is_DP", thisProperty)
133	+	PropertyMap(i)%is_DP = thisProperty
134	+
135	+	if (thisProperty) then
136	+	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
137	+	PropertyMap(i)%dipole_moment = thisDPproperty
138	+	endif
139	+
140	+	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
141	+	PropertyMap(i)%is_Sticky = thisProperty
142	+	call getElementProperty(atypes, i, "is_GB", thisProperty)
143	+	PropertyMap(i)%is_GB = thisProperty
144	+	call getElementProperty(atypes, i, "is_EAM", thisProperty)
145	+	PropertyMap(i)%is_EAM = thisProperty
146	+	end do
147	+
148	+	havePropertyMap = .true.
149	+
150	+	end subroutine createPropertyMap
151	+
152	+	subroutine setSimVariables()
153	+	SIM_uses_LJ = SimUsesLJ()
154	+	SIM_uses_sticky = SimUsesSticky()
155	+	SIM_uses_charges = SimUsesCharges()
156	+	SIM_uses_dipoles = SimUsesDipoles()
157	+	SIM_uses_RF = SimUsesRF()
158	+	SIM_uses_GB = SimUsesGB()
159	+	SIM_uses_EAM = SimUsesEAM()
160	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
161	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
162	+	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
163	+	SIM_uses_PBC = SimUsesPBC()
164	+	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
165	+
166	+	haveSIMvariables = .true.
167	+
168	+	return
169	+	end subroutine setSimVariables
170	+
171	+	subroutine doReadyCheck(error)
172	+	integer, intent(out) :: error
173	+
174	+	integer :: myStatus
175	+
176	+	error = 0
177	+
178	+	if (.not. havePropertyMap) then
179	+
180	+	myStatus = 0
181	+
182	+	call createPropertyMap(myStatus)
183	+
184	+	if (myStatus .ne. 0) then
185	+	write(default_error, *) 'createPropertyMap failed in do_Forces!'
186	+	error = -1
187	+	return
188	+	endif
189	+	endif
190	+
191	+	if (.not. haveSIMvariables) then
192	+	call setSimVariables()
193	+	endif
194	+
195	+	if (.not. haveRlist) then
196	+	write(default_error, *) 'rList has not been set in do_Forces!'
197	+	error = -1
198	+	return
199	+	endif
200	+
201	+	if (SIM_uses_LJ .and. FF_uses_LJ) then
202	+	if (.not. havePolicies) then
203	+	write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
204	+	error = -1
205	+	return
206	+	endif
207	+	endif
208
209	+	if (.not. haveNeighborList) then
210	+	write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
211	+	error = -1
212	+	return
213	+	end if
214	+
215	+	if (.not. haveSaneForceField) then
216	+	write(default_error, *) 'Force Field is not sane in do_Forces!'
217	+	error = -1
218	+	return
219	+	end if
220	+
221	+	#ifdef IS_MPI
222	+	if (.not. isMPISimSet()) then
223	+	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
224	+	error = -1
225	+	return
226	+	endif
227	+	#endif
228	+	return
229	+	end subroutine doReadyCheck
230	+
231	+
232		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
233
234		integer, intent(in) :: LJMIXPOLICY
#	Line 64 | Line 253 | contains
253
254		FF_uses_LJ = .false.
255		FF_uses_sticky = .false.
256	+	FF_uses_charges = .false.
257		FF_uses_dipoles = .false.
258		FF_uses_GB = .false.
259		FF_uses_EAM = .false.
260
261		call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
262		if (nMatches .gt. 0) FF_uses_LJ = .true.
263	<
263	>
264	>	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
265	>	if (nMatches .gt. 0) FF_uses_charges = .true.
266	>
267		call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
268		if (nMatches .gt. 0) FF_uses_dipoles = .true.
269
#	Line 84 | Line 277 | contains
277		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
278		if (nMatches .gt. 0) FF_uses_EAM = .true.
279
280	+	!! Assume sanity (for the sake of argument)
281	+	haveSaneForceField = .true.
282	+	!!
283	+	if (FF_uses_charges) then
284	+	dielect = getDielect()
285	+	call initialize_charge(dielect)
286	+	endif
287	+
288	+
289		!! check to make sure the FF_uses_RF setting makes sense
290
291		if (FF_uses_dipoles) then
90	–	rrf = getRrf()
91	–	rt = getRt()
92	–	call initialize_dipole(rrf, rt)
292		if (FF_uses_RF) then
293		dielect = getDielect()
294	<	call initialize_rf(rrf, rt, dielect)
294	>	call initialize_rf(dielect)
295		endif
296		else
297		if (FF_uses_RF) then
298		write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
299		thisStat = -1
300	+	haveSaneForceField = .false.
301		return
302		endif
303	<	endif
303	>	endif
304
305		if (FF_uses_LJ) then
306
107	–	call getRcut(rcut)
108	–
307		select case (LJMIXPOLICY)
308		case (LB_MIXING_RULE)
309	<	call init_lj_FF(LB_MIXING_RULE, rcut, my_status)
309	>	call init_lj_FF(LB_MIXING_RULE, my_status)
310		case (EXPLICIT_MIXING_RULE)
311	<	call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
311	>	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
312		case default
313		write(default_error,*) 'unknown LJ Mixing Policy!'
314		thisStat = -1
315	+	haveSaneForceField = .false.
316		return
317		end select
318		if (my_status /= 0) then
319		thisStat = -1
320	+	haveSaneForceField = .false.
321		return
322		end if
323	+	havePolicies = .true.
324		endif
325
326		if (FF_uses_sticky) then
327		call check_sticky_FF(my_status)
328		if (my_status /= 0) then
329		thisStat = -1
330	+	haveSaneForceField = .false.
331		return
332		end if
333		endif
334	<
334	>
335	>
336	>	if (FF_uses_EAM) then
337	>	call init_EAM_FF(my_status)
338	>	if (my_status /= 0) then
339	>	write(default_error, *) "init_EAM_FF returned a bad status"
340	>	thisStat = -1
341	>	haveSaneForceField = .false.
342	>	return
343	>	end if
344	>	endif
345	>
346		if (FF_uses_GB) then
347		call check_gb_pair_FF(my_status)
348		if (my_status .ne. 0) then
349		thisStat = -1
350	+	haveSaneForceField = .false.
351		return
352		endif
353		endif
354
355		if (FF_uses_GB .and. FF_uses_LJ) then
356		endif
357	<	if (.not. do_forces_initialized) then
357	>	if (.not. haveNeighborList) then
358		!! Create neighbor lists
359	<	call expandNeighborList(getNlocal(), my_status)
359	>	call expandNeighborList(nLocal, my_status)
360		if (my_Status /= 0) then
361		write(default_error,*) "SimSetup: ExpandNeighborList returned error."
362		thisStat = -1
363		return
364		endif
365	+	haveNeighborList = .true.
366		endif
367	<
153	<	do_forces_initialized = .true.
154	<
367	>
368		end subroutine init_FF
369
370
371		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
372		!------------------------------------------------------------->
373	<	subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
374	<	error)
373	>	subroutine do_force_loop(q, qcom, A, u_l, f, t, tau, pot, &
374	>	do_pot_c, do_stress_c, error)
375		!! Position array provided by C, dimensioned by getNlocal
376	<	real ( kind = dp ), dimension(3,getNlocal()) :: q
376	>	real ( kind = dp ), dimension(3,nLocal) :: q
377	>	!! molecular center-of-mass position array
378	>	real ( kind = dp ), dimension(3,nLocal) :: qcom
379		!! Rotation Matrix for each long range particle in simulation.
380	<	real( kind = dp), dimension(9,getNlocal()) :: A
380	>	real( kind = dp), dimension(9,nLocal) :: A
381		!! Unit vectors for dipoles (lab frame)
382	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
382	>	real( kind = dp ), dimension(3,nLocal) :: u_l
383		!! Force array provided by C, dimensioned by getNlocal
384	<	real ( kind = dp ), dimension(3,getNlocal()) :: f
384	>	real ( kind = dp ), dimension(3,nLocal) :: f
385		!! Torsion array provided by C, dimensioned by getNlocal
386	<	real( kind = dp ), dimension(3,getNlocal()) :: t
386	>	real( kind = dp ), dimension(3,nLocal) :: t
387	>
388		!! Stress Tensor
389		real( kind = dp), dimension(9) :: tau
390		real ( kind = dp ) :: pot
#	Line 179 | Line 395 | contains
395		real( kind = DP ) :: pot_local
396		integer :: nrow
397		integer :: ncol
398	+	integer :: nprocs
399		#endif
183	–	integer :: nlocal
400		integer :: natoms
401		logical :: update_nlist
402		integer :: i, j, jbeg, jend, jnab
403		integer :: nlist
404	<	real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
405	<	real(kind=dp),dimension(3) :: d
404	>	real( kind = DP ) ::  rijsq, rcijsq
405	>	real(kind=dp),dimension(3) :: d, dc
406		real(kind=dp) :: rfpot, mu_i, virial
407	<	integer :: me_i
407	>	integer :: me_i, me_j
408		logical :: is_dp_i
409		integer :: neighborListSize
410		integer :: listerror, error
411		integer :: localError
412	+	integer :: propPack_i, propPack_j
413
414	+	real(kind=dp) :: listSkin = 1.0
415	+
416		!! initialize local variables
417	<
417	>
418		#ifdef IS_MPI
419		pot_local = 0.0_dp
201	–	nlocal = getNlocal()
420		nrow   = getNrow(plan_row)
421		ncol   = getNcol(plan_col)
422		#else
205	–	nlocal = getNlocal()
423		natoms = nlocal
424		#endif
208	–
209	–	call getRcut(rcut,rc2=rcutsq)
210	–	call getRlist(rlist,rlistsq)
425
426	<	call check_initialization(localError)
426	>	call doReadyCheck(localError)
427		if ( localError .ne. 0 ) then
428	+	call handleError("do_force_loop", "Not Initialized")
429		error = -1
430		return
431		end if
432		call zero_work_arrays()
433	<
433	>
434		do_pot = do_pot_c
435		do_stress = do_stress_c
436
222	–
437		! Gather all information needed by all force loops:
438
439		#ifdef IS_MPI
440	<
440	>
441		call gather(q,q_Row,plan_row3d)
442		call gather(q,q_Col,plan_col3d)
443	<
444	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
443	>
444	>	if (SIM_uses_molecular_cutoffs) then
445	>	call gather(qcom,qcom_Row,plan_row3d)
446	>	call gather(qcom,qcom_Col,plan_col3d)
447	>	endif
448	>
449	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
450		call gather(u_l,u_l_Row,plan_row3d)
451		call gather(u_l,u_l_Col,plan_col3d)
452
#	Line 237 | Line 456 | contains
456
457		#endif
458
459	<	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
459	>	!! Begin force loop timing:
460	>	#ifdef PROFILE
461	>	call cpu_time(forceTimeInitial)
462	>	nloops = nloops + 1
463	>	#endif
464	>
465	>	if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
466		!! See if we need to update neighbor lists
467	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
467	>	if (SIM_uses_molecular_cutoffs) then
468	>	call checkNeighborList(nlocal, qcom, listSkin, update_nlist)
469	>	else
470	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
471	>	endif
472		!! if_mpi_gather_stuff_for_prepair
473		!! do_prepair_loop_if_needed
474		!! if_mpi_scatter_stuff_from_prepair
475		!! if_mpi_gather_stuff_from_prepair_to_main_loop
476	<	else
477	<	!! See if we need to update neighbor lists
478	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
476	>
477	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
478	>	#ifdef IS_MPI
479	>
480	>	if (update_nlist) then
481	>
482	>	!! save current configuration, construct neighbor list,
483	>	!! and calculate forces
484	>	if (SIM_uses_molecular_cutoffs) then
485	>	call saveNeighborList(nlocal, qcom)
486	>	else
487	>	call saveNeighborList(nlocal, q)
488	>	endif
489	>
490	>	neighborListSize = size(list)
491	>	nlist = 0
492	>
493	>	do i = 1, nrow
494	>	point(i) = nlist + 1
495	>
496	>	prepair_inner: do j = 1, ncol
497	>
498	>	if (skipThisPair(i,j)) cycle prepair_inner
499	>
500	>	if (SIM_uses_molecular_cutoffs) then
501	>	call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
502	>	dc, rcijsq)
503	>	else
504	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
505	>	dc, rcijsq)
506	>	endif
507	>
508	>	if (rcijsq < rlistsq) then
509	>
510	>	nlist = nlist + 1
511	>
512	>	if (nlist > neighborListSize) then
513	>	call expandNeighborList(nlocal, listerror)
514	>	if (listerror /= 0) then
515	>	error = -1
516	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
517	>	return
518	>	end if
519	>	neighborListSize = size(list)
520	>	endif
521	>
522	>	list(nlist) = j
523	>
524	>	if (SIM_uses_molecular_cutoffs) then
525	>	! since the simulation distances were in molecular COMs:
526	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
527	>	d, rijsq)
528	>	else
529	>	d(1:3) = dc(1:3)
530	>	rijsq = rcijsq
531	>	endif
532	>
533	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
534	>	do_pot, do_stress, u_l, A, f, t, pot_local)
535	>	endif
536	>	enddo prepair_inner
537	>	enddo
538	>
539	>	point(nrow + 1) = nlist + 1
540	>
541	>	else  !! (of update_check)
542	>
543	>	! use the list to find the neighbors
544	>	do i = 1, nrow
545	>	JBEG = POINT(i)
546	>	JEND = POINT(i+1) - 1
547	>	! check thiat molecule i has neighbors
548	>	if (jbeg .le. jend) then
549	>
550	>	do jnab = jbeg, jend
551	>	j = list(jnab)
552	>
553	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
554	>	d, rijsq)
555	>	if (SIM_uses_molecular_cutoffs) then
556	>	call get_interatomic_vector(qcom_Row(:,i),qcom_Col(:,j),&
557	>	dc, rcijsq)
558	>	else
559	>	dc(1:3) = d(1:3)
560	>	rcijsq = rijsq
561	>	endif
562	>
563	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
564	>	do_pot, do_stress, u_l, A, f, t, pot_local)
565	>
566	>	enddo
567	>	endif
568	>	enddo
569	>	endif
570	>
571	>	#else
572	>
573	>	if (update_nlist) then
574	>
575	>	! save current configuration, contruct neighbor list,
576	>	! and calculate forces
577	>	call saveNeighborList(natoms, q)
578	>
579	>	neighborListSize = size(list)
580	>
581	>	nlist = 0
582	>
583	>	do i = 1, natoms-1
584	>	point(i) = nlist + 1
585	>
586	>	prepair_inner: do j = i+1, natoms
587	>
588	>	if (skipThisPair(i,j))  cycle prepair_inner
589	>
590	>	if (SIM_uses_molecular_cutoffs) then
591	>	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
592	>	dc, rcijsq)
593	>	else
594	>	call get_interatomic_vector(q(:,i), q(:,j), dc, rcijsq)
595	>	endif
596	>
597	>	if (rcijsq < rlistsq) then
598	>
599	>	nlist = nlist + 1
600	>
601	>	if (nlist > neighborListSize) then
602	>	call expandNeighborList(natoms, listerror)
603	>	if (listerror /= 0) then
604	>	error = -1
605	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
606	>	return
607	>	end if
608	>	neighborListSize = size(list)
609	>	endif
610	>
611	>	list(nlist) = j
612	>
613	>
614	>	if (SIM_uses_molecular_cutoffs) then
615	>	! since the simulation distances were in molecular COMs:
616	>	call get_interatomic_vector(q(:,i), q(:,j), &
617	>	d, rijsq)
618	>	else
619	>	dc(1:3) = d(1:3)
620	>	rcijsq = rijsq
621	>	endif
622	>
623	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
624	>	do_pot, do_stress, u_l, A, f, t, pot)
625	>
626	>	endif
627	>	enddo prepair_inner
628	>	enddo
629	>
630	>	point(natoms) = nlist + 1
631	>
632	>	else !! (update)
633	>
634	>	! use the list to find the neighbors
635	>	do i = 1, natoms-1
636	>	JBEG = POINT(i)
637	>	JEND = POINT(i+1) - 1
638	>	! check thiat molecule i has neighbors
639	>	if (jbeg .le. jend) then
640	>
641	>	do jnab = jbeg, jend
642	>	j = list(jnab)
643	>
644	>	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
645	>	if (SIM_uses_molecular_cutoffs) then
646	>	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
647	>	dc, rcijsq)
648	>	else
649	>	dc(1:3) = d(1:3)
650	>	rcijsq = rijsq
651	>	endif
652	>
653	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
654	>	do_pot, do_stress, u_l, A, f, t, pot)
655	>
656	>	enddo
657	>	endif
658	>	enddo
659	>	endif
660	>	#endif
661	>	!! Do rest of preforce calculations
662	>	!! do necessary preforce calculations
663	>	call do_preforce(nlocal,pot)
664	>	! we have already updated the neighbor list set it to false...
665	>	update_nlist = .false.
666	>	else
667	>	!! See if we need to update neighbor lists for non pre-pair
668	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
669		endif
670
671	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
672	+
673		#ifdef IS_MPI
674
675		if (update_nlist) then
255	–
676		!! save current configuration, construct neighbor list,
677		!! and calculate forces
678		call saveNeighborList(nlocal, q)
#	Line 261 | Line 681 | contains
681		nlist = 0
682
683		do i = 1, nrow
684	+
685		point(i) = nlist + 1
686
687		inner: do j = 1, ncol
688
689		if (skipThisPair(i,j)) cycle inner
690
691	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
691	>	if (SIM_uses_molecular_cutoffs) then
692	>	call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
693	>	dc, rcijsq)
694	>	else
695	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
696	>	dc, rcijsq)
697	>	endif
698
699	<	if (rijsq <  rlistsq) then
699	>	if (rcijsq < rlistsq) then
700
701		nlist = nlist + 1
702
#	Line 284 | Line 711 | contains
711		endif
712
713		list(nlist) = j
714	<
715	<	if (rijsq <  rcutsq) then
716	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
717	<	u_l, A, f, t, pot_local)
714	>
715	>	if (SIM_uses_molecular_cutoffs) then
716	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
717	>	d, rijsq)
718	>	else
719	>	d(1:3) = dc(1:3)
720	>	rijsq = rcijsq
721		endif
722	+
723	+	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
724	+	do_pot, do_stress, u_l, A, f, t, pot_local)
725	+
726		endif
727		enddo inner
728		enddo
729	<
729	>
730		point(nrow + 1) = nlist + 1
731
732		else  !! (of update_check)
733	<
733	>
734		! use the list to find the neighbors
735		do i = 1, nrow
736		JBEG = POINT(i)
#	Line 306 | Line 740 | contains
740
741		do jnab = jbeg, jend
742		j = list(jnab)
743	<
743	>
744		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
745	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
746	<	u_l, A, f, t, pot_local)
747	<
745	>	if (SIM_uses_molecular_cutoffs) then
746	>	call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
747	>	dc, rcijsq)
748	>	else
749	>	dc(1:3) = d(1:3)
750	>	rcijsq = rijsq
751	>	endif
752	>
753	>	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
754	>	do_pot, do_stress, u_l, A, f, t, pot_local)
755	>
756		enddo
757		endif
758		enddo
#	Line 325 | Line 767 | contains
767		call saveNeighborList(natoms, q)
768
769		neighborListSize = size(list)
770	<
770	>
771		nlist = 0
772
773		do i = 1, natoms-1
#	Line 334 | Line 776 | contains
776		inner: do j = i+1, natoms
777
778		if (skipThisPair(i,j))  cycle inner
779	<
780	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
781	<
782	<
783	<	if (rijsq <  rlistsq) then
779	>
780	>	if (SIM_uses_molecular_cutoffs) then
781	>	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
782	>	dc, rcijsq)
783	>	else
784	>	call get_interatomic_vector(q(:,i), q(:,j), &
785	>	dc, rcijsq)
786	>	endif
787	>
788	>	if (rcijsq < rlistsq) then
789
790		nlist = nlist + 1
791	<
791	>
792		if (nlist > neighborListSize) then
793		call expandNeighborList(natoms, listerror)
794		if (listerror /= 0) then
#	Line 354 | Line 801 | contains
801
802		list(nlist) = j
803
804	<	if (rijsq <  rcutsq) then
805	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
806	<	u_l, A, f, t, pot)
804	>	if (SIM_uses_molecular_cutoffs) then
805	>	call get_interatomic_vector(q(:,i), q(:,j), &
806	>	d, rijsq)
807	>	else
808	>	d(1:3) = dc(1:3)
809	>	rijsq = rcijsq
810		endif
811	+
812	+	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
813	+	do_pot, do_stress, u_l, A, f, t, pot)
814	+
815		endif
816		enddo inner
817		enddo
#	Line 376 | Line 830 | contains
830		do jnab = jbeg, jend
831		j = list(jnab)
832
379	–	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
380	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
381	–	u_l, A, f, t, pot)
833
834	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
835	+	if (SIM_uses_molecular_cutoffs) then
836	+	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
837	+	dc, rcijsq)
838	+	else
839	+	dc(1:3) = d(1:3)
840	+	rcijsq = rijsq
841	+	endif
842	+
843	+	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
844	+	do_pot, do_stress, u_l, A, f, t, pot)
845	+
846		enddo
847		endif
848		enddo
#	Line 389 | Line 852 | contains
852
853		! phew, done with main loop.
854
855	+	!! Do timing
856	+	#ifdef PROFILE
857	+	call cpu_time(forceTimeFinal)
858	+	forceTime = forceTime + forceTimeFinal - forceTimeInitial
859	+	#endif
860	+
861	+
862		#ifdef IS_MPI
863		!!distribute forces
864	<
864	>
865		f_temp = 0.0_dp
866		call scatter(f_Row,f_temp,plan_row3d)
867		do i = 1,nlocal
868		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
869		end do
870	<
870	>
871		f_temp = 0.0_dp
872		call scatter(f_Col,f_temp,plan_col3d)
873		do i = 1,nlocal
874		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
875		end do
876
877	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
877	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
878		t_temp = 0.0_dp
879		call scatter(t_Row,t_temp,plan_row3d)
880		do i = 1,nlocal
#	Line 421 | Line 891 | contains
891		if (do_pot) then
892		! scatter/gather pot_row into the members of my column
893		call scatter(pot_Row, pot_Temp, plan_row)
894	<
894	>
895		! scatter/gather pot_local into all other procs
896		! add resultant to get total pot
897		do i = 1, nlocal
#	Line 429 | Line 899 | contains
899		enddo
900
901		pot_Temp = 0.0_DP
902	<
902	>
903		call scatter(pot_Col, pot_Temp, plan_col)
904		do i = 1, nlocal
905		pot_local = pot_local + pot_Temp(i)
906		enddo
907
908	<	endif
908	>	endif
909		#endif
910	<
911	<	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
910	>
911	>	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
912
913	<	if (FF_uses_RF .and. SimUsesRF()) then
913	>	if (FF_uses_RF .and. SIM_uses_RF) then
914
915		#ifdef IS_MPI
916		call scatter(rf_Row,rf,plan_row3d)
#	Line 450 | Line 920 | contains
920		end do
921		#endif
922
923	<	do i = 1, getNlocal()
924	<
923	>	do i = 1, nLocal
924	>
925		rfpot = 0.0_DP
926		#ifdef IS_MPI
927		me_i = atid_row(i)
928		#else
929		me_i = atid(i)
930		#endif
931	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
932	<	if ( is_DP_i ) then
933	<	call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
931	>
932	>	if (PropertyMap(me_i)%is_DP) then
933	>
934	>	mu_i = PropertyMap(me_i)%dipole_moment
935	>
936		!! The reaction field needs to include a self contribution
937		!! to the field:
938	<	call accumulate_self_rf(i, mu_i, u_l)
938	>	call accumulate_self_rf(i, mu_i, u_l)
939		!! Get the reaction field contribution to the
940		!! potential and torques:
941		call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#	Line 477 | Line 949 | contains
949		enddo
950		endif
951		endif
952	<
953	<
952	>
953	>
954		#ifdef IS_MPI
955	<
955	>
956		if (do_pot) then
957		pot = pot + pot_local
958		!! we assume the c code will do the allreduce to get the total potential
959		!! we could do it right here if we needed to...
960		endif
961	<
961	>
962		if (do_stress) then
963	<	call mpi_allreduce(tau_Temp, tau,9,mpi_double_precision,mpi_sum, &
963	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
964		mpi_comm_world,mpi_err)
965		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
966		mpi_comm_world,mpi_err)
967		endif
968	<
968	>
969		#else
970	<
970	>
971		if (do_stress) then
972		tau = tau_Temp
973		virial = virial_Temp
974		endif
975	<
975	>
976		#endif
977
978	+
979		end subroutine do_force_loop
980	+
981	+	subroutine do_pair(i, j, rijsq, d, rcijsq, dc, mfact, do_pot, do_stress, &
982	+	u_l, A, f, t, pot)
983
508	–	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
509	–
984		real( kind = dp ) :: pot
985	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
986	<	real (kind=dp), dimension(9,getNlocal()) :: A
987	<	real (kind=dp), dimension(3,getNlocal()) :: f
988	<	real (kind=dp), dimension(3,getNlocal()) :: t
985	>	real( kind = dp ), dimension(nLocal)   :: mfact
986	>	real( kind = dp ), dimension(3,nLocal) :: u_l
987	>	real( kind = dp ), dimension(9,nLocal) :: A
988	>	real( kind = dp ), dimension(3,nLocal) :: f
989	>	real( kind = dp ), dimension(3,nLocal) :: t
990
991		logical, intent(inout) :: do_pot, do_stress
992		integer, intent(in) :: i, j
993	<	real ( kind = dp ), intent(inout)    :: rijsq
994	<	real ( kind = dp )                :: r
995	<	real ( kind = dp ), intent(inout) :: d(3)
521	<	logical :: is_LJ_i, is_LJ_j
522	<	logical :: is_DP_i, is_DP_j
523	<	logical :: is_GB_i, is_GB_j
524	<	logical :: is_Sticky_i, is_Sticky_j
993	>	real ( kind = dp ), intent(inout) :: rijsq, rcijsq
994	>	real ( kind = dp )                :: r, rc
995	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
996		integer :: me_i, me_j
997
998		r = sqrt(rijsq)
999	+	if (SIM_uses_molecular_cutoffs) then
1000	+	rc = sqrt(rcijsq)
1001	+	else
1002	+	rc = r
1003	+	endif
1004
529	–	#ifdef IS_MPI
1005
1006	+	#ifdef IS_MPI
1007	+	if (tagRow(i) .eq. tagColumn(j)) then
1008	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
1009	+	endif
1010		me_i = atid_row(i)
1011		me_j = atid_col(j)
533	–
1012		#else
535	–
1013		me_i = atid(i)
1014		me_j = atid(j)
538	–
1015		#endif
1016	<
1017	<	if (FF_uses_LJ .and. SimUsesLJ()) then
1018	<	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
1019	<	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
1020	<
1021	<	if ( is_LJ_i .and. is_LJ_j ) &
1022	<	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1016	>
1017	>	if (FF_uses_LJ .and. SIM_uses_LJ) then
1018	>
1019	>	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
1020	>	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1021	>	endif
1022	>
1023		endif
1024	<
1025	<	if (FF_uses_dipoles .and. SimUsesDipoles()) then
550	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
551	<	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
1024	>
1025	>	if (FF_uses_charges .and. SIM_uses_charges) then
1026
1027	<	if ( is_DP_i .and. is_DP_j ) then
1028	<
1027	>	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
1028	>	call do_charge_pair(i, j, d, r, rijsq, dc, rc, rcijsq, mfact, &
1029	>	pot, f, do_pot, do_stress, SIM_uses_molecular_cutoffs)
1030	>	endif
1031	>
1032	>	endif
1033	>
1034	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
1035	>
1036	>	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
1037		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
1038		do_pot, do_stress)
1039	<	if (FF_uses_RF .and. SimUsesRF()) then
1039	>	if (FF_uses_RF .and. SIM_uses_RF) then
1040		call accumulate_rf(i, j, r, u_l)
1041		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
1042	<	endif
561	<
1042	>	endif
1043		endif
1044	+
1045		endif
1046
1047	<	if (FF_uses_Sticky .and. SimUsesSticky()) then
1047	>	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1048
1049	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
568	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
569	<
570	<	if ( is_Sticky_i .and. is_Sticky_j ) then
1049	>	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1050		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
1051		do_pot, do_stress)
1052		endif
1053	+
1054		endif
1055
1056
1057	<	if (FF_uses_GB .and. SimUsesGB()) then
578	<
579	<	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
580	<	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
1057	>	if (FF_uses_GB .and. SIM_uses_GB) then
1058
1059	<	if ( is_GB_i .and. is_GB_j ) then
1059	>	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1060		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
1061		do_pot, do_stress)
1062		endif
586	–	endif
587	–
588	–	end subroutine do_pair
1063
1064	<
591	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
592	<
593	<	real (kind = dp), dimension(3) :: q_i
594	<	real (kind = dp), dimension(3) :: q_j
595	<	real ( kind = dp ), intent(out) :: r_sq
596	<	real( kind = dp ) :: d(3)
597	<	real( kind = dp ) :: d_old(3)
598	<	d(1:3) = q_j(1:3) - q_i(1:3)
599	<	d_old = d
600	<	! Wrap back into periodic box if necessary
601	<	if ( SimUsesPBC() ) then
1064	>	endif
1065
1066	<	d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
604	<	int(abs(d(1:3)/box(1:3)) + 0.5_dp)
1066	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1067
1068	+	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1069	+	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1070	+	endif
1071	+
1072		endif
607	–	r_sq = dot_product(d,d)
608	–
609	–	end subroutine get_interatomic_vector
1073
1074	<	subroutine check_initialization(error)
1075	<	integer, intent(out) :: error
1076	<
1077	<	error = 0
1078	<	! Make sure we are properly initialized.
1079	<	if (.not. do_forces_initialized) then
1080	<	error = -1
1081	<	return
1074	>	end subroutine do_pair
1075	>
1076	>
1077	>
1078	>	subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, &
1079	>	do_pot, do_stress, u_l, A, f, t, pot)
1080	>	real( kind = dp ) :: pot
1081	>	real( kind = dp ), dimension(3,nLocal) :: u_l
1082	>	real (kind=dp), dimension(9,nLocal) :: A
1083	>	real (kind=dp), dimension(3,nLocal) :: f
1084	>	real (kind=dp), dimension(3,nLocal) :: t
1085	>
1086	>	logical, intent(inout) :: do_pot, do_stress
1087	>	integer, intent(in) :: i, j
1088	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1089	>	real ( kind = dp )                :: r, rc
1090	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1091	>
1092	>	logical :: is_EAM_i, is_EAM_j
1093	>
1094	>	integer :: me_i, me_j
1095	>
1096	>
1097	>	r = sqrt(rijsq)
1098	>	if (SIM_uses_molecular_cutoffs) then
1099	>	rc = sqrt(rcijsq)
1100	>	else
1101	>	rc = r
1102		endif
1103	+
1104
1105		#ifdef IS_MPI
1106	<	if (.not. isMPISimSet()) then
1107	<	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
1108	<	error = -1
1109	<	return
1110	<	endif
1106	>	if (tagRow(i) .eq. tagColumn(j)) then
1107	>	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
1108	>	endif
1109	>
1110	>	me_i = atid_row(i)
1111	>	me_j = atid_col(j)
1112	>
1113	>	#else
1114	>
1115	>	me_i = atid(i)
1116	>	me_j = atid(j)
1117	>
1118		#endif
1119
1120	<	return
630	<	end subroutine check_initialization
1120	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1121
1122	<
1123	<	subroutine zero_work_arrays()
634	<
635	<	#ifdef IS_MPI
1122	>	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1123	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1124
1125	<	q_Row = 0.0_dp
1126	<	q_Col = 0.0_dp
1127	<
640	<	u_l_Row = 0.0_dp
641	<	u_l_Col = 0.0_dp
642	<
643	<	A_Row = 0.0_dp
644	<	A_Col = 0.0_dp
645	<
646	<	f_Row = 0.0_dp
647	<	f_Col = 0.0_dp
648	<	f_Temp = 0.0_dp
649	<
650	<	t_Row = 0.0_dp
651	<	t_Col = 0.0_dp
652	<	t_Temp = 0.0_dp
653	<
654	<	pot_Row = 0.0_dp
655	<	pot_Col = 0.0_dp
656	<	pot_Temp = 0.0_dp
1125	>	endif
1126	>
1127	>	end subroutine do_prepair
1128
658	–	rf_Row = 0.0_dp
659	–	rf_Col = 0.0_dp
660	–	rf_Temp = 0.0_dp
1129
662	–	#endif
1130
664	–	rf = 0.0_dp
665	–	tau_Temp = 0.0_dp
666	–	virial_Temp = 0.0_dp
667	–	end subroutine zero_work_arrays
668	–
669	–	function skipThisPair(atom1, atom2) result(skip_it)
670	–	integer, intent(in) :: atom1
671	–	integer, intent(in), optional :: atom2
672	–	logical :: skip_it
673	–	integer :: unique_id_1, unique_id_2
674	–	integer :: me_i,me_j
675	–	integer :: i
1131
1132	<	skip_it = .false.
1133	<
1134	<	!! there are a number of reasons to skip a pair or a particle
1135	<	!! mostly we do this to exclude atoms who are involved in short
1136	<	!! range interactions (bonds, bends, torsions), but we also need
1137	<	!! to exclude some overcounted interactions that result from
1138	<	!! the parallel decomposition
1139	<
1140	<	#ifdef IS_MPI
1141	<	!! in MPI, we have to look up the unique IDs for each atom
1142	<	unique_id_1 = tagRow(atom1)
1143	<	#else
1144	<	!! in the normal loop, the atom numbers are unique
1145	<	unique_id_1 = atom1
1146	<	#endif
1147	<
1148	<	!! We were called with only one atom, so just check the global exclude
1149	<	!! list for this atom
1150	<	if (.not. present(atom2)) then
1151	<	do i = 1, nExcludes_global
1152	<	if (excludesGlobal(i) == unique_id_1) then
1153	<	skip_it = .true.
1154	<	return
1155	<	end if
1156	<	end do
1157	<	return
1158	<	end if
1159	<
1132	>	subroutine do_preforce(nlocal,pot)
1133	>	integer :: nlocal
1134	>	real( kind = dp ) :: pot
1135	>
1136	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1137	>	call calc_EAM_preforce_Frho(nlocal,pot)
1138	>	endif
1139	>
1140	>
1141	>	end subroutine do_preforce
1142	>
1143	>
1144	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1145	>
1146	>	real (kind = dp), dimension(3) :: q_i
1147	>	real (kind = dp), dimension(3) :: q_j
1148	>	real ( kind = dp ), intent(out) :: r_sq
1149	>	real( kind = dp ) :: d(3), scaled(3)
1150	>	integer i
1151	>
1152	>	d(1:3) = q_j(1:3) - q_i(1:3)
1153	>
1154	>	! Wrap back into periodic box if necessary
1155	>	if ( SIM_uses_PBC ) then
1156	>
1157	>	if( .not.boxIsOrthorhombic ) then
1158	>	! calc the scaled coordinates.
1159	>
1160	>	scaled = matmul(HmatInv, d)
1161	>
1162	>	! wrap the scaled coordinates
1163	>
1164	>	scaled = scaled  - anint(scaled)
1165	>
1166	>
1167	>	! calc the wrapped real coordinates from the wrapped scaled
1168	>	! coordinates
1169	>
1170	>	d = matmul(Hmat,scaled)
1171	>
1172	>	else
1173	>	! calc the scaled coordinates.
1174	>
1175	>	do i = 1, 3
1176	>	scaled(i) = d(i) * HmatInv(i,i)
1177	>
1178	>	! wrap the scaled coordinates
1179	>
1180	>	scaled(i) = scaled(i) - anint(scaled(i))
1181	>
1182	>	! calc the wrapped real coordinates from the wrapped scaled
1183	>	! coordinates
1184	>
1185	>	d(i) = scaled(i)*Hmat(i,i)
1186	>	enddo
1187	>	endif
1188	>
1189	>	endif
1190	>
1191	>	r_sq = dot_product(d,d)
1192	>
1193	>	end subroutine get_interatomic_vector
1194	>
1195	>	subroutine zero_work_arrays()
1196	>
1197		#ifdef IS_MPI
1198	<	unique_id_2 = tagColumn(atom2)
1198	>
1199	>	q_Row = 0.0_dp
1200	>	q_Col = 0.0_dp
1201	>
1202	>	qcom_Row = 0.0_dp
1203	>	qcom_Col = 0.0_dp
1204	>
1205	>	u_l_Row = 0.0_dp
1206	>	u_l_Col = 0.0_dp
1207	>
1208	>	A_Row = 0.0_dp
1209	>	A_Col = 0.0_dp
1210	>
1211	>	f_Row = 0.0_dp
1212	>	f_Col = 0.0_dp
1213	>	f_Temp = 0.0_dp
1214	>
1215	>	t_Row = 0.0_dp
1216	>	t_Col = 0.0_dp
1217	>	t_Temp = 0.0_dp
1218	>
1219	>	pot_Row = 0.0_dp
1220	>	pot_Col = 0.0_dp
1221	>	pot_Temp = 0.0_dp
1222	>
1223	>	rf_Row = 0.0_dp
1224	>	rf_Col = 0.0_dp
1225	>	rf_Temp = 0.0_dp
1226	>
1227	>	#endif
1228	>
1229	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1230	>	call clean_EAM()
1231	>	endif
1232	>
1233	>	rf = 0.0_dp
1234	>	tau_Temp = 0.0_dp
1235	>	virial_Temp = 0.0_dp
1236	>	end subroutine zero_work_arrays
1237	>
1238	>	function skipThisPair(atom1, atom2) result(skip_it)
1239	>	integer, intent(in) :: atom1
1240	>	integer, intent(in), optional :: atom2
1241	>	logical :: skip_it
1242	>	integer :: unique_id_1, unique_id_2
1243	>	integer :: me_i,me_j
1244	>	integer :: i
1245	>
1246	>	skip_it = .false.
1247	>
1248	>	!! there are a number of reasons to skip a pair or a particle
1249	>	!! mostly we do this to exclude atoms who are involved in short
1250	>	!! range interactions (bonds, bends, torsions), but we also need
1251	>	!! to exclude some overcounted interactions that result from
1252	>	!! the parallel decomposition
1253	>
1254	>	#ifdef IS_MPI
1255	>	!! in MPI, we have to look up the unique IDs for each atom
1256	>	unique_id_1 = tagRow(atom1)
1257		#else
1258	<	unique_id_2 = atom2
1258	>	!! in the normal loop, the atom numbers are unique
1259	>	unique_id_1 = atom1
1260		#endif
1261	<
1261	>
1262	>	!! We were called with only one atom, so just check the global exclude
1263	>	!! list for this atom
1264	>	if (.not. present(atom2)) then
1265	>	do i = 1, nExcludes_global
1266	>	if (excludesGlobal(i) == unique_id_1) then
1267	>	skip_it = .true.
1268	>	return
1269	>	end if
1270	>	end do
1271	>	return
1272	>	end if
1273	>
1274		#ifdef IS_MPI
1275	<	!! this situation should only arise in MPI simulations
1276	<	if (unique_id_1 == unique_id_2) then
1277	<	skip_it = .true.
715	<	return
716	<	end if
717	<
718	<	!! this prevents us from doing the pair on multiple processors
719	<	if (unique_id_1 < unique_id_2) then
720	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
721	<	skip_it = .true.
722	<	return
723	<	endif
724	<	else
725	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
726	<	skip_it = .true.
727	<	return
728	<	endif
729	<	endif
1275	>	unique_id_2 = tagColumn(atom2)
1276	>	#else
1277	>	unique_id_2 = atom2
1278		#endif
1279	+
1280	+	#ifdef IS_MPI
1281	+	!! this situation should only arise in MPI simulations
1282	+	if (unique_id_1 == unique_id_2) then
1283	+	skip_it = .true.
1284	+	return
1285	+	end if
1286	+
1287	+	!! this prevents us from doing the pair on multiple processors
1288	+	if (unique_id_1 < unique_id_2) then
1289	+	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1290	+	skip_it = .true.
1291	+	return
1292	+	endif
1293	+	else
1294	+	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1295	+	skip_it = .true.
1296	+	return
1297	+	endif
1298	+	endif
1299	+	#endif
1300	+
1301	+	!! the rest of these situations can happen in all simulations:
1302	+	do i = 1, nExcludes_global
1303	+	if ((excludesGlobal(i) == unique_id_1) .or. &
1304	+	(excludesGlobal(i) == unique_id_2)) then
1305	+	skip_it = .true.
1306	+	return
1307	+	endif
1308	+	enddo
1309	+
1310	+	do i = 1, nExcludes_local
1311	+	if (excludesLocal(1,i) == unique_id_1) then
1312	+	if (excludesLocal(2,i) == unique_id_2) then
1313	+	skip_it = .true.
1314	+	return
1315	+	endif
1316	+	else
1317	+	if (excludesLocal(1,i) == unique_id_2) then
1318	+	if (excludesLocal(2,i) == unique_id_1) then
1319	+	skip_it = .true.
1320	+	return
1321	+	endif
1322	+	endif
1323	+	endif
1324	+	end do
1325	+
1326	+	return
1327	+	end function skipThisPair
1328
1329	<	!! the rest of these situations can happen in all simulations:
1330	<	do i = 1, nExcludes_global
1331	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1332	<	(excludesGlobal(i) == unique_id_2)) then
1333	<	skip_it = .true.
1334	<	return
1335	<	endif
1336	<	enddo
1337	<
1338	<	do i = 1, nExcludes_local
1339	<	if (excludesLocal(1,i) == unique_id_1) then
1340	<	if (excludesLocal(2,i) == unique_id_2) then
1341	<	skip_it = .true.
1342	<	return
1343	<	endif
1344	<	else
1345	<	if (excludesLocal(1,i) == unique_id_2) then
1346	<	if (excludesLocal(2,i) == unique_id_1) then
1347	<	skip_it = .true.
1348	<	return
1349	<	endif
1350	<	endif
1351	<	endif
1352	<	end do
1353	<
757	<	return
758	<	end function skipThisPair
759	<
760	<	function FF_UsesDirectionalAtoms() result(doesit)
761	<	logical :: doesit
762	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
763	<	FF_uses_GB .or. FF_uses_RF
764	<	end function FF_UsesDirectionalAtoms
765	<
766	<	function FF_RequiresPrepairCalc() result(doesit)
767	<	logical :: doesit
768	<	doesit = FF_uses_EAM
769	<	end function FF_RequiresPrepairCalc
770	<
771	<	function FF_RequiresPostpairCalc() result(doesit)
772	<	logical :: doesit
773	<	doesit = FF_uses_RF
774	<	end function FF_RequiresPostpairCalc
775	<
1329	>	function FF_UsesDirectionalAtoms() result(doesit)
1330	>	logical :: doesit
1331	>	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1332	>	FF_uses_GB .or. FF_uses_RF
1333	>	end function FF_UsesDirectionalAtoms
1334	>
1335	>	function FF_RequiresPrepairCalc() result(doesit)
1336	>	logical :: doesit
1337	>	doesit = FF_uses_EAM
1338	>	end function FF_RequiresPrepairCalc
1339	>
1340	>	function FF_RequiresPostpairCalc() result(doesit)
1341	>	logical :: doesit
1342	>	doesit = FF_uses_RF
1343	>	end function FF_RequiresPostpairCalc
1344	>
1345	>	#ifdef PROFILE
1346	>	function getforcetime() result(totalforcetime)
1347	>	real(kind=dp) :: totalforcetime
1348	>	totalforcetime = forcetime
1349	>	end function getforcetime
1350	>	#endif
1351	>
1352	>	!! This cleans componets of force arrays belonging only to fortran
1353	>
1354		end module do_Forces

Diff Legend

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