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 572 by mmeineke, Wed Jul 2 21:26:55 2003 UTC vs.
Revision 1138 by gezelter, Wed Apr 28 21:39:22 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.17 2003-07-02 21:26:55 mmeineke Exp $, $Date: 2003-07-02 21:26:55 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $
7	>	!! @version $Id: do_Forces.F90,v 1.52 2004-04-28 21:39:22 gezelter Exp $, $Date: 2004-04-28 21:39:22 $, $Name: not supported by cvs2svn $, $Revision: 1.52 $
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, mfact, 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	>	!! mass factors used for molecular cutoffs
380	>	real ( kind = dp ), dimension(3,nLocal) :: mfact
381		!! Rotation Matrix for each long range particle in simulation.
382	<	real( kind = dp), dimension(9,getNlocal()) :: A
382	>	real( kind = dp), dimension(9,nLocal) :: A
383		!! Unit vectors for dipoles (lab frame)
384	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
384	>	real( kind = dp ), dimension(3,nLocal) :: u_l
385		!! Force array provided by C, dimensioned by getNlocal
386	<	real ( kind = dp ), dimension(3,getNlocal()) :: f
386	>	real ( kind = dp ), dimension(3,nLocal) :: f
387		!! Torsion array provided by C, dimensioned by getNlocal
388	<	real( kind = dp ), dimension(3,getNlocal()) :: t
388	>	real( kind = dp ), dimension(3,nLocal) :: t
389	>
390		!! Stress Tensor
391		real( kind = dp), dimension(9) :: tau
392		real ( kind = dp ) :: pot
#	Line 179 | Line 397 | contains
397		real( kind = DP ) :: pot_local
398		integer :: nrow
399		integer :: ncol
400	+	integer :: nprocs
401		#endif
183	–	integer :: nlocal
402		integer :: natoms
403		logical :: update_nlist
404		integer :: i, j, jbeg, jend, jnab
405		integer :: nlist
406	<	real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
407	<	real(kind=dp),dimension(3) :: d
406	>	real( kind = DP ) ::  rijsq, rcijsq
407	>	real(kind=dp),dimension(3) :: d, dc
408		real(kind=dp) :: rfpot, mu_i, virial
409	<	integer :: me_i
409	>	integer :: me_i, me_j
410		logical :: is_dp_i
411		integer :: neighborListSize
412		integer :: listerror, error
413		integer :: localError
414	+	integer :: propPack_i, propPack_j
415
416	+	real(kind=dp) :: listSkin = 1.0
417	+
418		!! initialize local variables
419	<
419	>
420		#ifdef IS_MPI
421		pot_local = 0.0_dp
201	–	nlocal = getNlocal()
422		nrow   = getNrow(plan_row)
423		ncol   = getNcol(plan_col)
424		#else
205	–	nlocal = getNlocal()
425		natoms = nlocal
426		#endif
208	–
209	–	call getRcut(rcut,rc2=rcutsq)
210	–	call getRlist(rlist,rlistsq)
427
428	<	call check_initialization(localError)
428	>	call doReadyCheck(localError)
429		if ( localError .ne. 0 ) then
430	+	call handleError("do_force_loop", "Not Initialized")
431		error = -1
432		return
433		end if
434		call zero_work_arrays()
435	<
435	>
436		do_pot = do_pot_c
437		do_stress = do_stress_c
438	<
438	>
439		! Gather all information needed by all force loops:
440
441		#ifdef IS_MPI
442	<
442	>
443		call gather(q,q_Row,plan_row3d)
444		call gather(q,q_Col,plan_col3d)
445	<
446	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
445	>
446	>	if (SIM_uses_molecular_cutoffs) then
447	>	call gather(qcom,qcom_Row,plan_row3d)
448	>	call gather(qcom,qcom_Col,plan_col3d)
449	>	endif
450	>
451	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
452		call gather(u_l,u_l_Row,plan_row3d)
453		call gather(u_l,u_l_Col,plan_col3d)
454
#	Line 236 | Line 458 | contains
458
459		#endif
460
461	<	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
461	>	!! Begin force loop timing:
462	>	#ifdef PROFILE
463	>	call cpu_time(forceTimeInitial)
464	>	nloops = nloops + 1
465	>	#endif
466	>
467	>	if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
468		!! See if we need to update neighbor lists
469	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
469	>	if (SIM_uses_molecular_cutoffs) then
470	>	call checkNeighborList(nlocal, qcom, listSkin, update_nlist)
471	>	else
472	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
473	>	endif
474		!! if_mpi_gather_stuff_for_prepair
475		!! do_prepair_loop_if_needed
476		!! if_mpi_scatter_stuff_from_prepair
477		!! if_mpi_gather_stuff_from_prepair_to_main_loop
478	<	else
479	<	!! See if we need to update neighbor lists
480	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
478	>
479	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
480	>	#ifdef IS_MPI
481	>
482	>	if (update_nlist) then
483	>
484	>	!! save current configuration, construct neighbor list,
485	>	!! and calculate forces
486	>	if (SIM_uses_molecular_cutoffs) then
487	>	call saveNeighborList(nlocal, qcom)
488	>	else
489	>	call saveNeighborList(nlocal, q)
490	>	endif
491	>
492	>	neighborListSize = size(list)
493	>	nlist = 0
494	>
495	>	do i = 1, nrow
496	>	point(i) = nlist + 1
497	>
498	>	prepair_inner: do j = 1, ncol
499	>
500	>	if (skipThisPair(i,j)) cycle prepair_inner
501	>
502	>	if (SIM_uses_molecular_cutoffs) then
503	>	call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
504	>	dc, rcijsq)
505	>	else
506	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
507	>	dc, rcijsq)
508	>	endif
509	>
510	>	if (rcijsq < rlistsq) then
511	>
512	>	nlist = nlist + 1
513	>
514	>	if (nlist > neighborListSize) then
515	>	call expandNeighborList(nlocal, listerror)
516	>	if (listerror /= 0) then
517	>	error = -1
518	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
519	>	return
520	>	end if
521	>	neighborListSize = size(list)
522	>	endif
523	>
524	>	list(nlist) = j
525	>
526	>	if (SIM_uses_molecular_cutoffs) then
527	>	! since the simulation distances were in molecular COMs:
528	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
529	>	d, rijsq)
530	>	else
531	>	d(1:3) = dc(1:3)
532	>	rijsq = rcijsq
533	>	endif
534	>
535	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
536	>	do_pot, do_stress, u_l, A, f, t, pot_local)
537	>	endif
538	>	enddo prepair_inner
539	>	enddo
540	>
541	>	point(nrow + 1) = nlist + 1
542	>
543	>	else  !! (of update_check)
544	>
545	>	! use the list to find the neighbors
546	>	do i = 1, nrow
547	>	JBEG = POINT(i)
548	>	JEND = POINT(i+1) - 1
549	>	! check thiat molecule i has neighbors
550	>	if (jbeg .le. jend) then
551	>
552	>	do jnab = jbeg, jend
553	>	j = list(jnab)
554	>
555	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
556	>	d, rijsq)
557	>	if (SIM_uses_molecular_cutoffs) then
558	>	call get_interatomic_vector(qcom_Row(:,i),qcom_Col(:,j),&
559	>	dc, rcijsq)
560	>	else
561	>	dc(1:3) = d(1:3)
562	>	rcijsq = rijsq
563	>	endif
564	>
565	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
566	>	do_pot, do_stress, u_l, A, f, t, pot_local)
567	>
568	>	enddo
569	>	endif
570	>	enddo
571	>	endif
572	>
573	>	#else
574	>
575	>	if (update_nlist) then
576	>
577	>	! save current configuration, contruct neighbor list,
578	>	! and calculate forces
579	>	call saveNeighborList(natoms, q)
580	>
581	>	neighborListSize = size(list)
582	>
583	>	nlist = 0
584	>
585	>	do i = 1, natoms-1
586	>	point(i) = nlist + 1
587	>
588	>	prepair_inner: do j = i+1, natoms
589	>
590	>	if (skipThisPair(i,j))  cycle prepair_inner
591	>
592	>	if (SIM_uses_molecular_cutoffs) then
593	>	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
594	>	dc, rcijsq)
595	>	else
596	>	call get_interatomic_vector(q(:,i), q(:,j), dc, rcijsq)
597	>	endif
598	>
599	>	if (rcijsq < rlistsq) then
600	>
601	>	nlist = nlist + 1
602	>
603	>	if (nlist > neighborListSize) then
604	>	call expandNeighborList(natoms, listerror)
605	>	if (listerror /= 0) then
606	>	error = -1
607	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
608	>	return
609	>	end if
610	>	neighborListSize = size(list)
611	>	endif
612	>
613	>	list(nlist) = j
614	>
615	>
616	>	if (SIM_uses_molecular_cutoffs) then
617	>	! since the simulation distances were in molecular COMs:
618	>	call get_interatomic_vector(q(:,i), q(:,j), &
619	>	d, rijsq)
620	>	else
621	>	dc(1:3) = d(1:3)
622	>	rcijsq = rijsq
623	>	endif
624	>
625	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
626	>	do_pot, do_stress, u_l, A, f, t, pot)
627	>
628	>	endif
629	>	enddo prepair_inner
630	>	enddo
631	>
632	>	point(natoms) = nlist + 1
633	>
634	>	else !! (update)
635	>
636	>	! use the list to find the neighbors
637	>	do i = 1, natoms-1
638	>	JBEG = POINT(i)
639	>	JEND = POINT(i+1) - 1
640	>	! check thiat molecule i has neighbors
641	>	if (jbeg .le. jend) then
642	>
643	>	do jnab = jbeg, jend
644	>	j = list(jnab)
645	>
646	>	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
647	>	if (SIM_uses_molecular_cutoffs) then
648	>	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
649	>	dc, rcijsq)
650	>	else
651	>	dc(1:3) = d(1:3)
652	>	rcijsq = rijsq
653	>	endif
654	>
655	>	call do_prepair(i, j, rijsq, d, rcijsq, dc, &
656	>	do_pot, do_stress, u_l, A, f, t, pot)
657	>
658	>	enddo
659	>	endif
660	>	enddo
661	>	endif
662	>	#endif
663	>	!! Do rest of preforce calculations
664	>	!! do necessary preforce calculations
665	>	call do_preforce(nlocal,pot)
666	>	! we have already updated the neighbor list set it to false...
667	>	update_nlist = .false.
668	>	else
669	>	!! See if we need to update neighbor lists for non pre-pair
670	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
671		endif
672
673	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
674	+
675		#ifdef IS_MPI
676
677		if (update_nlist) then
254	–
678		!! save current configuration, construct neighbor list,
679		!! and calculate forces
680		call saveNeighborList(nlocal, q)
#	Line 260 | Line 683 | contains
683		nlist = 0
684
685		do i = 1, nrow
686	+
687		point(i) = nlist + 1
688
689		inner: do j = 1, ncol
690
691		if (skipThisPair(i,j)) cycle inner
692
693	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
693	>	if (SIM_uses_molecular_cutoffs) then
694	>	call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
695	>	dc, rcijsq)
696	>	else
697	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
698	>	dc, rcijsq)
699	>	endif
700
701	<	if (rijsq <  rlistsq) then
701	>	if (rcijsq < rlistsq) then
702
703		nlist = nlist + 1
704
#	Line 283 | Line 713 | contains
713		endif
714
715		list(nlist) = j
716	<
717	<	if (rijsq <  rcutsq) then
718	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
719	<	u_l, A, f, t, pot_local)
716	>
717	>	if (SIM_uses_molecular_cutoffs) then
718	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
719	>	d, rijsq)
720	>	else
721	>	d(1:3) = dc(1:3)
722	>	rijsq = rcijsq
723		endif
724	+
725	+	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
726	+	do_pot, do_stress, u_l, A, f, t, pot_local)
727	+
728		endif
729		enddo inner
730		enddo
731	<
731	>
732		point(nrow + 1) = nlist + 1
733
734		else  !! (of update_check)
735	<
735	>
736		! use the list to find the neighbors
737		do i = 1, nrow
738		JBEG = POINT(i)
#	Line 305 | Line 742 | contains
742
743		do jnab = jbeg, jend
744		j = list(jnab)
745	<
745	>
746		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
747	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
748	<	u_l, A, f, t, pot_local)
749	<
747	>	if (SIM_uses_molecular_cutoffs) then
748	>	call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
749	>	dc, rcijsq)
750	>	else
751	>	dc(1:3) = d(1:3)
752	>	rcijsq = rijsq
753	>	endif
754	>
755	>	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
756	>	do_pot, do_stress, u_l, A, f, t, pot_local)
757	>
758		enddo
759		endif
760		enddo
#	Line 324 | Line 769 | contains
769		call saveNeighborList(natoms, q)
770
771		neighborListSize = size(list)
772	<
772	>
773		nlist = 0
774
775		do i = 1, natoms-1
#	Line 333 | Line 778 | contains
778		inner: do j = i+1, natoms
779
780		if (skipThisPair(i,j))  cycle inner
781	<
782	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
783	<
784	<
785	<	if (rijsq <  rlistsq) then
781	>
782	>	if (SIM_uses_molecular_cutoffs) then
783	>	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
784	>	dc, rcijsq)
785	>	else
786	>	call get_interatomic_vector(q(:,i), q(:,j), &
787	>	dc, rcijsq)
788	>	endif
789	>
790	>	if (rcijsq < rlistsq) then
791
792		nlist = nlist + 1
793	<
793	>
794		if (nlist > neighborListSize) then
795		call expandNeighborList(natoms, listerror)
796		if (listerror /= 0) then
#	Line 353 | Line 803 | contains
803
804		list(nlist) = j
805
806	<	if (rijsq <  rcutsq) then
807	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
808	<	u_l, A, f, t, pot)
806	>	if (SIM_uses_molecular_cutoffs) then
807	>	call get_interatomic_vector(q(:,i), q(:,j), &
808	>	d, rijsq)
809	>	else
810	>	d(1:3) = dc(1:3)
811	>	rijsq = rcijsq
812		endif
813	+
814	+	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
815	+	do_pot, do_stress, u_l, A, f, t, pot)
816	+
817		endif
818		enddo inner
819		enddo
#	Line 375 | Line 832 | contains
832		do jnab = jbeg, jend
833		j = list(jnab)
834
378	–	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
379	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
380	–	u_l, A, f, t, pot)
835
836	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
837	+	if (SIM_uses_molecular_cutoffs) then
838	+	call get_interatomic_vector(qcom(:,i), qcom(:,j), &
839	+	dc, rcijsq)
840	+	else
841	+	dc(1:3) = d(1:3)
842	+	rcijsq = rijsq
843	+	endif
844	+
845	+	call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
846	+	do_pot, do_stress, u_l, A, f, t, pot)
847	+
848		enddo
849		endif
850		enddo
#	Line 388 | Line 854 | contains
854
855		! phew, done with main loop.
856
857	+	!! Do timing
858	+	#ifdef PROFILE
859	+	call cpu_time(forceTimeFinal)
860	+	forceTime = forceTime + forceTimeFinal - forceTimeInitial
861	+	#endif
862	+
863	+
864		#ifdef IS_MPI
865		!!distribute forces
866	<
866	>
867		f_temp = 0.0_dp
868		call scatter(f_Row,f_temp,plan_row3d)
869		do i = 1,nlocal
870		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
871		end do
872	<
872	>
873		f_temp = 0.0_dp
874		call scatter(f_Col,f_temp,plan_col3d)
875		do i = 1,nlocal
876		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
877		end do
878
879	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
879	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
880		t_temp = 0.0_dp
881		call scatter(t_Row,t_temp,plan_row3d)
882		do i = 1,nlocal
#	Line 420 | Line 893 | contains
893		if (do_pot) then
894		! scatter/gather pot_row into the members of my column
895		call scatter(pot_Row, pot_Temp, plan_row)
896	<
896	>
897		! scatter/gather pot_local into all other procs
898		! add resultant to get total pot
899		do i = 1, nlocal
#	Line 428 | Line 901 | contains
901		enddo
902
903		pot_Temp = 0.0_DP
904	<
904	>
905		call scatter(pot_Col, pot_Temp, plan_col)
906		do i = 1, nlocal
907		pot_local = pot_local + pot_Temp(i)
908		enddo
909
910	<	endif
910	>	endif
911		#endif
912	<
913	<	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
912	>
913	>	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
914
915	<	if (FF_uses_RF .and. SimUsesRF()) then
915	>	if (FF_uses_RF .and. SIM_uses_RF) then
916
917		#ifdef IS_MPI
918		call scatter(rf_Row,rf,plan_row3d)
#	Line 449 | Line 922 | contains
922		end do
923		#endif
924
925	<	do i = 1, getNlocal()
926	<
925	>	do i = 1, nLocal
926	>
927		rfpot = 0.0_DP
928		#ifdef IS_MPI
929		me_i = atid_row(i)
930		#else
931		me_i = atid(i)
932		#endif
933	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
934	<	if ( is_DP_i ) then
935	<	call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
933	>
934	>	if (PropertyMap(me_i)%is_DP) then
935	>
936	>	mu_i = PropertyMap(me_i)%dipole_moment
937	>
938		!! The reaction field needs to include a self contribution
939		!! to the field:
940	<	call accumulate_self_rf(i, mu_i, u_l)
940	>	call accumulate_self_rf(i, mu_i, u_l)
941		!! Get the reaction field contribution to the
942		!! potential and torques:
943		call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#	Line 476 | Line 951 | contains
951		enddo
952		endif
953		endif
954	<
955	<
954	>
955	>
956		#ifdef IS_MPI
957	<
957	>
958		if (do_pot) then
959		pot = pot + pot_local
960		!! we assume the c code will do the allreduce to get the total potential
961		!! we could do it right here if we needed to...
962		endif
963	<
963	>
964		if (do_stress) then
965	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
965	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
966		mpi_comm_world,mpi_err)
967		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
968		mpi_comm_world,mpi_err)
969		endif
970	<
970	>
971		#else
972	<
972	>
973		if (do_stress) then
974		tau = tau_Temp
975		virial = virial_Temp
976		endif
977	<
977	>
978		#endif
979
980	+
981		end subroutine do_force_loop
982	+
983	+	subroutine do_pair(i, j, rijsq, d, rcijsq, dc, mfact, do_pot, do_stress, &
984	+	u_l, A, f, t, pot)
985
507	–	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
508	–
986		real( kind = dp ) :: pot
987	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
988	<	real (kind=dp), dimension(9,getNlocal()) :: A
989	<	real (kind=dp), dimension(3,getNlocal()) :: f
990	<	real (kind=dp), dimension(3,getNlocal()) :: t
987	>	real( kind = dp ), dimension(nLocal)   :: mfact
988	>	real( kind = dp ), dimension(3,nLocal) :: u_l
989	>	real( kind = dp ), dimension(9,nLocal) :: A
990	>	real( kind = dp ), dimension(3,nLocal) :: f
991	>	real( kind = dp ), dimension(3,nLocal) :: t
992
993		logical, intent(inout) :: do_pot, do_stress
994		integer, intent(in) :: i, j
995	<	real ( kind = dp ), intent(inout)    :: rijsq
996	<	real ( kind = dp )                :: r
997	<	real ( kind = dp ), intent(inout) :: d(3)
520	<	logical :: is_LJ_i, is_LJ_j
521	<	logical :: is_DP_i, is_DP_j
522	<	logical :: is_GB_i, is_GB_j
523	<	logical :: is_Sticky_i, is_Sticky_j
995	>	real ( kind = dp ), intent(inout) :: rijsq, rcijsq
996	>	real ( kind = dp )                :: r, rc
997	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
998		integer :: me_i, me_j
999
1000		r = sqrt(rijsq)
1001	+	if (SIM_uses_molecular_cutoffs) then
1002	+	rc = sqrt(rcijsq)
1003	+	else
1004	+	rc = r
1005	+	endif
1006
1007
529	–
1008		#ifdef IS_MPI
1009		if (tagRow(i) .eq. tagColumn(j)) then
1010		write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
1011		endif
534	–
1012		me_i = atid_row(i)
1013		me_j = atid_col(j)
537	–
1014		#else
539	–
1015		me_i = atid(i)
1016		me_j = atid(j)
542	–
1017		#endif
1018	<
1019	<	if (FF_uses_LJ .and. SimUsesLJ()) then
1020	<	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
1021	<	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
1022	<
1023	<	if ( is_LJ_i .and. is_LJ_j ) &
1024	<	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1018	>
1019	>	if (FF_uses_LJ .and. SIM_uses_LJ) then
1020	>
1021	>	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
1022	>	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1023	>	endif
1024	>
1025		endif
1026	<
1027	<	if (FF_uses_dipoles .and. SimUsesDipoles()) then
554	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
555	<	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
1026	>
1027	>	if (FF_uses_charges .and. SIM_uses_charges) then
1028
1029	<	if ( is_DP_i .and. is_DP_j ) then
1030	<
1029	>	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
1030	>	call do_charge_pair(i, j, d, r, rijsq, dc, rc, rcijsq, mfact, &
1031	>	pot, f, do_pot, do_stress, SIM_uses_molecular_cutoffs)
1032	>	endif
1033	>
1034	>	endif
1035	>
1036	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
1037	>
1038	>	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
1039		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
1040		do_pot, do_stress)
1041	<	if (FF_uses_RF .and. SimUsesRF()) then
1041	>	if (FF_uses_RF .and. SIM_uses_RF) then
1042		call accumulate_rf(i, j, r, u_l)
1043		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
1044	<	endif
565	<
1044	>	endif
1045		endif
1046	+
1047		endif
1048
1049	<	if (FF_uses_Sticky .and. SimUsesSticky()) then
1049	>	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1050
1051	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
572	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
573	<
574	<	if ( is_Sticky_i .and. is_Sticky_j ) then
1051	>	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1052		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
1053		do_pot, do_stress)
1054		endif
1055	+
1056		endif
1057
1058
1059	<	if (FF_uses_GB .and. SimUsesGB()) then
582	<
583	<	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
584	<	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
1059	>	if (FF_uses_GB .and. SIM_uses_GB) then
1060
1061	<	if ( is_GB_i .and. is_GB_j ) then
1061	>	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1062		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
1063		do_pot, do_stress)
1064		endif
590	–	endif
591	–
592	–	end subroutine do_pair
1065
1066	<
595	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
596	<
597	<	real (kind = dp), dimension(3) :: q_i
598	<	real (kind = dp), dimension(3) :: q_j
599	<	real ( kind = dp ), intent(out) :: r_sq
600	<	real( kind = dp ) :: d(3), scaled(3)
601	<	integer i
602	<
603	<	d(1:3) = q_j(1:3) - q_i(1:3)
604	<
605	<	! Wrap back into periodic box if necessary
606	<	if ( SimUsesPBC() ) then
1066	>	endif
1067
1068	<	if( .not.boxIsOrthorhombic ) then
1069	<	! calc the scaled coordinates.
1070	<
1071	<	scaled = matmul(HmatInv, d)
1072	<
1073	<	! wrap the scaled coordinates
1068	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1069	>
1070	>	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1071	>	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1072	>	endif
1073	>
1074	>	endif
1075
1076	<	scaled = scaled  - anint(scaled)
616	<
1076	>	end subroutine do_pair
1077
618	–	! calc the wrapped real coordinates from the wrapped scaled
619	–	! coordinates
1078
621	–	d = matmul(Hmat,scaled)
1079
1080	<	else
1081	<	! calc the scaled coordinates.
1082	<
1083	<	do i = 1, 3
1084	<	scaled(i) = d(i) * HmatInv(i,i)
1085	<
1086	<	! wrap the scaled coordinates
1087	<
1088	<	scaled(i) = scaled(i) - anint(scaled(i))
1089	<
1090	<	! calc the wrapped real coordinates from the wrapped scaled
1091	<	! coordinates
1080	>	subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, &
1081	>	do_pot, do_stress, u_l, A, f, t, pot)
1082	>	real( kind = dp ) :: pot
1083	>	real( kind = dp ), dimension(3,nLocal) :: u_l
1084	>	real (kind=dp), dimension(9,nLocal) :: A
1085	>	real (kind=dp), dimension(3,nLocal) :: f
1086	>	real (kind=dp), dimension(3,nLocal) :: t
1087	>
1088	>	logical, intent(inout) :: do_pot, do_stress
1089	>	integer, intent(in) :: i, j
1090	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1091	>	real ( kind = dp )                :: r, rc
1092	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1093	>
1094	>	logical :: is_EAM_i, is_EAM_j
1095	>
1096	>	integer :: me_i, me_j
1097	>
1098
1099	<	d(i) = scaled(i)*Hmat(i,i)
1100	<	enddo
1101	<	endif
1102	<
1099	>	r = sqrt(rijsq)
1100	>	if (SIM_uses_molecular_cutoffs) then
1101	>	rc = sqrt(rcijsq)
1102	>	else
1103	>	rc = r
1104		endif
1105	<
642	<	r_sq = dot_product(d,d)
643	<
644	<	end subroutine get_interatomic_vector
645	<
646	<	subroutine check_initialization(error)
647	<	integer, intent(out) :: error
648	<
649	<	error = 0
650	<	! Make sure we are properly initialized.
651	<	if (.not. do_forces_initialized) then
652	<	error = -1
653	<	return
654	<	endif
1105	>
1106
1107		#ifdef IS_MPI
1108	<	if (.not. isMPISimSet()) then
1109	<	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
1110	<	error = -1
1111	<	return
1112	<	endif
1108	>	if (tagRow(i) .eq. tagColumn(j)) then
1109	>	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
1110	>	endif
1111	>
1112	>	me_i = atid_row(i)
1113	>	me_j = atid_col(j)
1114	>
1115	>	#else
1116	>
1117	>	me_i = atid(i)
1118	>	me_j = atid(j)
1119	>
1120		#endif
1121
1122	<	return
665	<	end subroutine check_initialization
1122	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1123
1124	<
1125	<	subroutine zero_work_arrays()
669	<
670	<	#ifdef IS_MPI
1124	>	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1125	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1126
1127	<	q_Row = 0.0_dp
1128	<	q_Col = 0.0_dp
1129	<
675	<	u_l_Row = 0.0_dp
676	<	u_l_Col = 0.0_dp
677	<
678	<	A_Row = 0.0_dp
679	<	A_Col = 0.0_dp
680	<
681	<	f_Row = 0.0_dp
682	<	f_Col = 0.0_dp
683	<	f_Temp = 0.0_dp
684	<
685	<	t_Row = 0.0_dp
686	<	t_Col = 0.0_dp
687	<	t_Temp = 0.0_dp
688	<
689	<	pot_Row = 0.0_dp
690	<	pot_Col = 0.0_dp
691	<	pot_Temp = 0.0_dp
1127	>	endif
1128	>
1129	>	end subroutine do_prepair
1130
693	–	rf_Row = 0.0_dp
694	–	rf_Col = 0.0_dp
695	–	rf_Temp = 0.0_dp
1131
697	–	#endif
1132
699	–	rf = 0.0_dp
700	–	tau_Temp = 0.0_dp
701	–	virial_Temp = 0.0_dp
702	–	end subroutine zero_work_arrays
703	–
704	–	function skipThisPair(atom1, atom2) result(skip_it)
705	–	integer, intent(in) :: atom1
706	–	integer, intent(in), optional :: atom2
707	–	logical :: skip_it
708	–	integer :: unique_id_1, unique_id_2
709	–	integer :: me_i,me_j
710	–	integer :: i
1133
1134	<	skip_it = .false.
1135	<
1136	<	!! there are a number of reasons to skip a pair or a particle
1137	<	!! mostly we do this to exclude atoms who are involved in short
1138	<	!! range interactions (bonds, bends, torsions), but we also need
1139	<	!! to exclude some overcounted interactions that result from
1140	<	!! the parallel decomposition
1141	<
1142	<	#ifdef IS_MPI
1143	<	!! in MPI, we have to look up the unique IDs for each atom
1144	<	unique_id_1 = tagRow(atom1)
1145	<	#else
1146	<	!! in the normal loop, the atom numbers are unique
1147	<	unique_id_1 = atom1
1148	<	#endif
1134	>	subroutine do_preforce(nlocal,pot)
1135	>	integer :: nlocal
1136	>	real( kind = dp ) :: pot
1137	>
1138	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1139	>	call calc_EAM_preforce_Frho(nlocal,pot)
1140	>	endif
1141	>
1142	>
1143	>	end subroutine do_preforce
1144	>
1145	>
1146	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1147	>
1148	>	real (kind = dp), dimension(3) :: q_i
1149	>	real (kind = dp), dimension(3) :: q_j
1150	>	real ( kind = dp ), intent(out) :: r_sq
1151	>	real( kind = dp ) :: d(3), scaled(3)
1152	>	integer i
1153	>
1154	>	d(1:3) = q_j(1:3) - q_i(1:3)
1155	>
1156	>	! Wrap back into periodic box if necessary
1157	>	if ( SIM_uses_PBC ) then
1158	>
1159	>	if( .not.boxIsOrthorhombic ) then
1160	>	! calc the scaled coordinates.
1161	>
1162	>	scaled = matmul(HmatInv, d)
1163	>
1164	>	! wrap the scaled coordinates
1165	>
1166	>	scaled = scaled  - anint(scaled)
1167	>
1168	>
1169	>	! calc the wrapped real coordinates from the wrapped scaled
1170	>	! coordinates
1171	>
1172	>	d = matmul(Hmat,scaled)
1173	>
1174	>	else
1175	>	! calc the scaled coordinates.
1176	>
1177	>	do i = 1, 3
1178	>	scaled(i) = d(i) * HmatInv(i,i)
1179	>
1180	>	! wrap the scaled coordinates
1181	>
1182	>	scaled(i) = scaled(i) - anint(scaled(i))
1183	>
1184	>	! calc the wrapped real coordinates from the wrapped scaled
1185	>	! coordinates
1186	>
1187	>	d(i) = scaled(i)*Hmat(i,i)
1188	>	enddo
1189	>	endif
1190	>
1191	>	endif
1192	>
1193	>	r_sq = dot_product(d,d)
1194	>
1195	>	end subroutine get_interatomic_vector
1196	>
1197	>	subroutine zero_work_arrays()
1198	>
1199	>	#ifdef IS_MPI
1200	>
1201	>	q_Row = 0.0_dp
1202	>	q_Col = 0.0_dp
1203
1204	<	!! We were called with only one atom, so just check the global exclude
1205	<	!! list for this atom
1206	<	if (.not. present(atom2)) then
1207	<	do i = 1, nExcludes_global
1208	<	if (excludesGlobal(i) == unique_id_1) then
1209	<	skip_it = .true.
1210	<	return
1211	<	end if
1212	<	end do
1213	<	return
1214	<	end if
1215	<
1204	>	qcom_Row = 0.0_dp
1205	>	qcom_Col = 0.0_dp
1206	>
1207	>	u_l_Row = 0.0_dp
1208	>	u_l_Col = 0.0_dp
1209	>
1210	>	A_Row = 0.0_dp
1211	>	A_Col = 0.0_dp
1212	>
1213	>	f_Row = 0.0_dp
1214	>	f_Col = 0.0_dp
1215	>	f_Temp = 0.0_dp
1216	>
1217	>	t_Row = 0.0_dp
1218	>	t_Col = 0.0_dp
1219	>	t_Temp = 0.0_dp
1220	>
1221	>	pot_Row = 0.0_dp
1222	>	pot_Col = 0.0_dp
1223	>	pot_Temp = 0.0_dp
1224	>
1225	>	rf_Row = 0.0_dp
1226	>	rf_Col = 0.0_dp
1227	>	rf_Temp = 0.0_dp
1228	>
1229	>	#endif
1230	>
1231	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1232	>	call clean_EAM()
1233	>	endif
1234	>
1235	>	rf = 0.0_dp
1236	>	tau_Temp = 0.0_dp
1237	>	virial_Temp = 0.0_dp
1238	>	end subroutine zero_work_arrays
1239	>
1240	>	function skipThisPair(atom1, atom2) result(skip_it)
1241	>	integer, intent(in) :: atom1
1242	>	integer, intent(in), optional :: atom2
1243	>	logical :: skip_it
1244	>	integer :: unique_id_1, unique_id_2
1245	>	integer :: me_i,me_j
1246	>	integer :: i
1247	>
1248	>	skip_it = .false.
1249	>
1250	>	!! there are a number of reasons to skip a pair or a particle
1251	>	!! mostly we do this to exclude atoms who are involved in short
1252	>	!! range interactions (bonds, bends, torsions), but we also need
1253	>	!! to exclude some overcounted interactions that result from
1254	>	!! the parallel decomposition
1255	>
1256		#ifdef IS_MPI
1257	<	unique_id_2 = tagColumn(atom2)
1257	>	!! in MPI, we have to look up the unique IDs for each atom
1258	>	unique_id_1 = tagRow(atom1)
1259		#else
1260	<	unique_id_2 = atom2
1260	>	!! in the normal loop, the atom numbers are unique
1261	>	unique_id_1 = atom1
1262		#endif
1263	<
1263	>
1264	>	!! We were called with only one atom, so just check the global exclude
1265	>	!! list for this atom
1266	>	if (.not. present(atom2)) then
1267	>	do i = 1, nExcludes_global
1268	>	if (excludesGlobal(i) == unique_id_1) then
1269	>	skip_it = .true.
1270	>	return
1271	>	end if
1272	>	end do
1273	>	return
1274	>	end if
1275	>
1276		#ifdef IS_MPI
1277	<	!! this situation should only arise in MPI simulations
1278	<	if (unique_id_1 == unique_id_2) then
1279	<	skip_it = .true.
750	<	return
751	<	end if
752	<
753	<	!! this prevents us from doing the pair on multiple processors
754	<	if (unique_id_1 < unique_id_2) then
755	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
756	<	skip_it = .true.
757	<	return
758	<	endif
759	<	else
760	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
761	<	skip_it = .true.
762	<	return
763	<	endif
764	<	endif
1277	>	unique_id_2 = tagColumn(atom2)
1278	>	#else
1279	>	unique_id_2 = atom2
1280		#endif
1281	+
1282	+	#ifdef IS_MPI
1283	+	!! this situation should only arise in MPI simulations
1284	+	if (unique_id_1 == unique_id_2) then
1285	+	skip_it = .true.
1286	+	return
1287	+	end if
1288	+
1289	+	!! this prevents us from doing the pair on multiple processors
1290	+	if (unique_id_1 < unique_id_2) then
1291	+	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1292	+	skip_it = .true.
1293	+	return
1294	+	endif
1295	+	else
1296	+	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1297	+	skip_it = .true.
1298	+	return
1299	+	endif
1300	+	endif
1301	+	#endif
1302	+
1303	+	!! the rest of these situations can happen in all simulations:
1304	+	do i = 1, nExcludes_global
1305	+	if ((excludesGlobal(i) == unique_id_1) .or. &
1306	+	(excludesGlobal(i) == unique_id_2)) then
1307	+	skip_it = .true.
1308	+	return
1309	+	endif
1310	+	enddo
1311	+
1312	+	do i = 1, nExcludes_local
1313	+	if (excludesLocal(1,i) == unique_id_1) then
1314	+	if (excludesLocal(2,i) == unique_id_2) then
1315	+	skip_it = .true.
1316	+	return
1317	+	endif
1318	+	else
1319	+	if (excludesLocal(1,i) == unique_id_2) then
1320	+	if (excludesLocal(2,i) == unique_id_1) then
1321	+	skip_it = .true.
1322	+	return
1323	+	endif
1324	+	endif
1325	+	endif
1326	+	end do
1327	+
1328	+	return
1329	+	end function skipThisPair
1330
1331	<	!! the rest of these situations can happen in all simulations:
1332	<	do i = 1, nExcludes_global
1333	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1334	<	(excludesGlobal(i) == unique_id_2)) then
1335	<	skip_it = .true.
1336	<	return
1337	<	endif
1338	<	enddo
1339	<
1340	<	do i = 1, nExcludes_local
1341	<	if (excludesLocal(1,i) == unique_id_1) then
1342	<	if (excludesLocal(2,i) == unique_id_2) then
1343	<	skip_it = .true.
1344	<	return
1345	<	endif
1346	<	else
1347	<	if (excludesLocal(1,i) == unique_id_2) then
1348	<	if (excludesLocal(2,i) == unique_id_1) then
1349	<	skip_it = .true.
1350	<	return
1351	<	endif
1352	<	endif
1353	<	endif
1354	<	end do
1355	<
792	<	return
793	<	end function skipThisPair
794	<
795	<	function FF_UsesDirectionalAtoms() result(doesit)
796	<	logical :: doesit
797	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
798	<	FF_uses_GB .or. FF_uses_RF
799	<	end function FF_UsesDirectionalAtoms
800	<
801	<	function FF_RequiresPrepairCalc() result(doesit)
802	<	logical :: doesit
803	<	doesit = FF_uses_EAM
804	<	end function FF_RequiresPrepairCalc
805	<
806	<	function FF_RequiresPostpairCalc() result(doesit)
807	<	logical :: doesit
808	<	doesit = FF_uses_RF
809	<	end function FF_RequiresPostpairCalc
810	<
1331	>	function FF_UsesDirectionalAtoms() result(doesit)
1332	>	logical :: doesit
1333	>	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1334	>	FF_uses_GB .or. FF_uses_RF
1335	>	end function FF_UsesDirectionalAtoms
1336	>
1337	>	function FF_RequiresPrepairCalc() result(doesit)
1338	>	logical :: doesit
1339	>	doesit = FF_uses_EAM
1340	>	end function FF_RequiresPrepairCalc
1341	>
1342	>	function FF_RequiresPostpairCalc() result(doesit)
1343	>	logical :: doesit
1344	>	doesit = FF_uses_RF
1345	>	end function FF_RequiresPostpairCalc
1346	>
1347	>	#ifdef PROFILE
1348	>	function getforcetime() result(totalforcetime)
1349	>	real(kind=dp) :: totalforcetime
1350	>	totalforcetime = forcetime
1351	>	end function getforcetime
1352	>	#endif
1353	>
1354	>	!! This cleans componets of force arrays belonging only to fortran
1355	>
1356		end module do_Forces

Diff Legend

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