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 843 by mmeineke, Wed Oct 29 20:41:39 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.35 2003-10-29 20:41:39 mmeineke Exp $, $Date: 2003-10-29 20:41:39 $, $Name: not supported by cvs2svn $, $Revision: 1.35 $
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
#	Line 30 | Line 31 | module do_Forces
31		#define __FORTRAN90
32		#include "fForceField.h"
33
34	<	logical, save :: do_forces_initialized = .false., haveRlist = .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
#	Line 45 | Line 63 | module do_Forces
63		public :: do_force_loop
64		public :: setRlistDF
65
66	+
67		#ifdef PROFILE
68	<	real(kind = dp) :: forceTime
69	<	real(kind = dp) :: forceTimeInitial, forceTimeFinal
70	<	real(kind = dp) :: globalForceTime
71	<	real(kind = dp) :: maxForceTime
53	<	integer, save :: nloops = 0
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 )
#	Line 63 | Line 94 | contains
94		rlistsq = rlist * rlist
95
96		haveRlist = .true.
66	–	if( havePolicies ) do_forces_initialized = .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 91 | 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 111 | 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
#	Line 122 | Line 297 | contains
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
#	Line 136 | Line 312 | contains
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
#	Line 156 | Line 336 | contains
336		if (FF_uses_EAM) then
337		call init_EAM_FF(my_status)
338		if (my_status /= 0) then
339	<	write(*,*) "init_EAM_FF returned a bad status"
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
165	–
166	–
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
187	–
188	–	havePolicies = .true.
189	–	if( haveRlist ) do_forces_initialized = .true.
190	–
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 217 | Line 399 | contains
399		integer :: ncol
400		integer :: nprocs
401		#endif
220	–	integer :: nlocal
402		integer :: natoms
403		logical :: update_nlist
404		integer :: i, j, jbeg, jend, jnab
405		integer :: nlist
406	<	real( kind = DP ) ::  rijsq
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
416	>	real(kind=dp) :: listSkin = 1.0
417
236	–
418		!! initialize local variables
419	<
419	>
420		#ifdef IS_MPI
421		pot_local = 0.0_dp
241	–	nlocal = getNlocal()
422		nrow   = getNrow(plan_row)
423		ncol   = getNcol(plan_col)
424		#else
245	–	nlocal = getNlocal()
425		natoms = nlocal
426		#endif
427	<
428	<	call check_initialization(localError)
427	>
428	>	call doReadyCheck(localError)
429		if ( localError .ne. 0 ) then
430	<	call handleError("do_force_loop","Not Initialized")
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	<
260	<
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 274 | Line 457 | contains
457		endif
458
459		#endif
460	<
461	<	!! Begin force loop timing:
460	>
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. SimRequiresPrepairCalc()) then
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, listSkin, 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	<
479	<	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
478	>
479	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
480		#ifdef IS_MPI
294	–
295	–	if (update_nlist) then
481
482	<	!! save current configuration, construct neighbor list,
298	<	!! and calculate forces
299	<	call saveNeighborList(nlocal, q)
300	<
301	<	neighborListSize = size(list)
302	<	nlist = 0
303	<
304	<	do i = 1, nrow
305	<	point(i) = nlist + 1
482	>	if (update_nlist) then
483
484	<	prepair_inner: do j = 1, ncol
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	<	if (skipThisPair(i,j)) cycle prepair_inner
310	<
311	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
312	<
313	<	if (rijsq < rlistsq) then
498	>	prepair_inner: do j = 1, ncol
499
500	<	nlist = nlist + 1
500	>	if (skipThisPair(i,j)) cycle prepair_inner
501
502	<	if (nlist > neighborListSize) then
503	<	call expandNeighborList(nlocal, listerror)
504	<	if (listerror /= 0) then
505	<	error = -1
506	<	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
507	<	return
323	<	end if
324	<	neighborListSize = size(list)
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	<	list(nlist) = j
511	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
512	<	endif
513	<	enddo prepair_inner
514	<	enddo
515	<
516	<	point(nrow + 1) = nlist + 1
517	<
518	<	else  !! (of update_check)
519	<
520	<	! use the list to find the neighbors
521	<	do i = 1, nrow
522	<	JBEG = POINT(i)
523	<	JEND = POINT(i+1) - 1
524	<	! check thiat molecule i has neighbors
525	<	if (jbeg .le. jend) then
526	<
527	<	do jnab = jbeg, jend
528	<	j = list(jnab)
529	<
530	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
531	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
532	<	u_l, A, f, t, pot_local)
533	<
534	<	enddo
535	<	endif
536	<	enddo
537	<	endif
538	<
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
357	–
358	–	if (update_nlist) then
574
575	<	! save current configuration, contruct neighbor list,
361	<	! and calculate forces
362	<	call saveNeighborList(natoms, q)
363	<
364	<	neighborListSize = size(list)
365	<
366	<	nlist = 0
367	<
368	<	do i = 1, natoms-1
369	<	point(i) = nlist + 1
575	>	if (update_nlist) then
576
577	<	prepair_inner: do j = i+1, natoms
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	<	if (skipThisPair(i,j))  cycle prepair_inner
374	<
375	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
376	<
377	<
378	<	if (rijsq < rlistsq) then
379	<
380	<
381	<	nlist = nlist + 1
382	<
383	<	if (nlist > neighborListSize) then
384	<	call expandNeighborList(natoms, listerror)
385	<	if (listerror /= 0) then
386	<	error = -1
387	<	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
388	<	return
389	<	end if
390	<	neighborListSize = size(list)
391	<	endif
588	>	prepair_inner: do j = i+1, natoms
589
590	<	list(nlist) = j
590	>	if (skipThisPair(i,j))  cycle prepair_inner
591
592	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
593	<	u_l, A, f, t, pot)
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	<	endif
600	<	enddo prepair_inner
601	<	enddo
602	<
603	<	point(natoms) = nlist + 1
604	<
605	<	else !! (update)
606	<
607	<	! use the list to find the neighbors
608	<	do i = 1, natoms-1
609	<	JBEG = POINT(i)
610	<	JEND = POINT(i+1) - 1
611	<	! check thiat molecule i has neighbors
612	<	if (jbeg .le. jend) then
613	<
614	<	do jnab = jbeg, jend
615	<	j = list(jnab)
616	<
617	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
618	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
619	<	u_l, A, f, t, pot)
620	<
621	<	enddo
622	<	endif
623	<	enddo
624	<	endif
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.
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	<
674	<
437	<
438	<
439	<	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
440	<
441	<
442	<
443	<
444	<
672	>
673	>	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
674	>
675		#ifdef IS_MPI
676
677		if (update_nlist) then
448	–
678		!! save current configuration, construct neighbor list,
679		!! and calculate forces
680		call saveNeighborList(nlocal, q)
#	Line 454 | 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 477 | Line 713 | contains
713		endif
714
715		list(nlist) = j
480	–
481	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
482	–	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 498 | 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 517 | 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 526 | 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 546 | Line 803 | contains
803
804		list(nlist) = j
805
806	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
807	<	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 567 | Line 832 | contains
832		do jnab = jbeg, jend
833		j = list(jnab)
834
570	–	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
571	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
572	–	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 579 | Line 853 | contains
853		#endif
854
855		! phew, done with main loop.
856	<
857	<	!! Do timing
856	>
857	>	!! Do timing
858		#ifdef PROFILE
859		call cpu_time(forceTimeFinal)
860		forceTime = forceTime + forceTimeFinal - forceTimeInitial
861		#endif
862	<
863	<
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 619 | 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 627 | 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 648 | 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 675 | 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	<	#ifdef PROFILE
705	<	if (do_pot) then
706	<
707	<	#ifdef IS_MPI
708	<
709	<
710	<	call printCommTime()
711	<
712	<	call mpi_allreduce(forceTime,globalForceTime,1,MPI_DOUBLE_PRECISION, &
713	<	mpi_sum,mpi_comm_world,mpi_err)
714	<
715	<	call mpi_allreduce(forceTime,maxForceTime,1,MPI_DOUBLE_PRECISION, &
716	<	MPI_MAX,mpi_comm_world,mpi_err)
717	<
718	<	call mpi_comm_size( MPI_COMM_WORLD, nprocs,mpi_err)
719	<
720	<	if (getMyNode() == 0) then
721	<	write(*,*) "Total processor time spent in force calculations is: ", globalForceTime
722	<	write(*,*) "Total Time spent in force loop per processor is: ", globalforceTime/nprocs
723	<	write(*,*) "Maximum force time on any processor is: ", maxForceTime
724	<	end if
725	<	#else
726	<	write(*,*) "Time spent in force loop is: ", forceTime
727	<	#endif
728	<
729	<
730	<	endif
731	<
732	<	#endif
733	<
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
736	–	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
737	–
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)
749	<	logical :: is_LJ_i, is_LJ_j
750	<	logical :: is_DP_i, is_DP_j
751	<	logical :: is_GB_i, is_GB_j
752	<	logical :: is_EAM_i,is_EAM_j
753	<	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	+
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
762	–
1012		me_i = atid_row(i)
1013		me_j = atid_col(j)
765	–
1014		#else
767	–
1015		me_i = atid(i)
1016		me_j = atid(j)
770	–
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
782	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
783	<	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
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
792	<
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)
799	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
800	<
801	<	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
809	<
810	<
811	<	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
812	<	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
818	–	endif
819	–
1065
1066	<
1067	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1068	<	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
1069	<	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
1070	<
1071	<	if ( is_EAM_i .and. is_EAM_j ) &
1072	<	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1073	<	endif
1074	<
1066	>	endif
1067	>
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
831	–
832	–
1076		end subroutine do_pair
1077
1078
1079
1080	<	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
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,getNlocal()) :: u_l
1084	<	real (kind=dp), dimension(9,getNlocal()) :: A
1085	<	real (kind=dp), dimension(3,getNlocal()) :: f
1086	<	real (kind=dp), dimension(3,getNlocal()) :: t
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
1091	<	real ( kind = dp )                :: r
1092	<	real ( kind = dp ), intent(inout) :: d(3)
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	<	r = sqrt(rijsq)
1098	>
1099	>	r = sqrt(rijsq)
1100	>	if (SIM_uses_molecular_cutoffs) then
1101	>	rc = sqrt(rcijsq)
1102	>	else
1103	>	rc = r
1104	>	endif
1105
1106
1107		#ifdef IS_MPI
#	Line 869 | Line 1119 | contains
1119
1120		#endif
1121
1122	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1123	<	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
1124	<	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
875	<
876	<	if ( is_EAM_i .and. is_EAM_j ) &
1122	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1123	>
1124	>	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1125		call calc_EAM_prepair_rho(i, j, d, r, rijsq )
878	–	endif
1126
1127	+	endif
1128	+
1129		end subroutine do_prepair
1130
1131
1132
1133
1134	<	subroutine do_preforce(nlocal,pot)
1135	<	integer :: nlocal
1136	<	real( kind = dp ) :: pot
1137	<
1138	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1139	<	call calc_EAM_preforce_Frho(nlocal,pot)
1140	<	endif
1141	<
1142	<
1143	<	end subroutine do_preforce
895	<
896	<
897	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
898	<
899	<	real (kind = dp), dimension(3) :: q_i
900	<	real (kind = dp), dimension(3) :: q_j
901	<	real ( kind = dp ), intent(out) :: r_sq
902	<	real( kind = dp ) :: d(3), scaled(3)
903	<	integer i
904	<
905	<	d(1:3) = q_j(1:3) - q_i(1:3)
906	<
907	<	! Wrap back into periodic box if necessary
908	<	if ( SimUsesPBC() ) then
909	<
910	<	if( .not.boxIsOrthorhombic ) then
911	<	! calc the scaled coordinates.
912	<
913	<	scaled = matmul(HmatInv, d)
914	<
915	<	! wrap the scaled coordinates
916	<
917	<	scaled = scaled  - anint(scaled)
918	<
919	<
920	<	! calc the wrapped real coordinates from the wrapped scaled
921	<	! coordinates
922	<
923	<	d = matmul(Hmat,scaled)
924	<
925	<	else
926	<	! calc the scaled coordinates.
927	<
928	<	do i = 1, 3
929	<	scaled(i) = d(i) * HmatInv(i,i)
930	<
931	<	! wrap the scaled coordinates
932	<
933	<	scaled(i) = scaled(i) - anint(scaled(i))
934	<
935	<	! calc the wrapped real coordinates from the wrapped scaled
936	<	! coordinates
937	<
938	<	d(i) = scaled(i)*Hmat(i,i)
939	<	enddo
940	<	endif
941	<
942	<	endif
943	<
944	<	r_sq = dot_product(d,d)
945	<
946	<	end subroutine get_interatomic_vector
947	<
948	<	subroutine check_initialization(error)
949	<	integer, intent(out) :: error
950	<
951	<	error = 0
952	<	! Make sure we are properly initialized.
953	<	if (.not. do_forces_initialized) then
954	<	write(*,*) "Forces not initialized"
955	<	error = -1
956	<	return
957	<	endif
958	<
959	<	#ifdef IS_MPI
960	<	if (.not. isMPISimSet()) then
961	<	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
962	<	error = -1
963	<	return
964	<	endif
965	<	#endif
966	<
967	<	return
968	<	end subroutine check_initialization
969	<
970	<
971	<	subroutine zero_work_arrays()
972	<
973	<	#ifdef IS_MPI
974	<
975	<	q_Row = 0.0_dp
976	<	q_Col = 0.0_dp
977	<
978	<	u_l_Row = 0.0_dp
979	<	u_l_Col = 0.0_dp
980	<
981	<	A_Row = 0.0_dp
982	<	A_Col = 0.0_dp
983	<
984	<	f_Row = 0.0_dp
985	<	f_Col = 0.0_dp
986	<	f_Temp = 0.0_dp
987	<
988	<	t_Row = 0.0_dp
989	<	t_Col = 0.0_dp
990	<	t_Temp = 0.0_dp
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	<	pot_Row = 0.0_dp
1146	<	pot_Col = 0.0_dp
1147	<	pot_Temp = 0.0_dp
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	<	rf_Row = 0.0_dp
1205	<	rf_Col = 0.0_dp
1206	<	rf_Temp = 0.0_dp
1207	<
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
1001	–
1230
1231	<	if (FF_uses_EAM .and. SimUsesEAM()) then
1232	<	call clean_EAM()
1233	<	endif
1234	<
1235	<
1236	<
1237	<
1238	<
1239	<	rf = 0.0_dp
1240	<	tau_Temp = 0.0_dp
1241	<	virial_Temp = 0.0_dp
1242	<	end subroutine zero_work_arrays
1243	<
1244	<	function skipThisPair(atom1, atom2) result(skip_it)
1245	<	integer, intent(in) :: atom1
1246	<	integer, intent(in), optional :: atom2
1247	<	logical :: skip_it
1248	<	integer :: unique_id_1, unique_id_2
1249	<	integer :: me_i,me_j
1250	<	integer :: i
1251	<
1252	<	skip_it = .false.
1253	<
1254	<	!! there are a number of reasons to skip a pair or a particle
1255	<	!! mostly we do this to exclude atoms who are involved in short
1028	<	!! range interactions (bonds, bends, torsions), but we also need
1029	<	!! to exclude some overcounted interactions that result from
1030	<	!! the parallel decomposition
1031	<
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	<	!! in MPI, we have to look up the unique IDs for each atom
1258	<	unique_id_1 = tagRow(atom1)
1257	>	!! in MPI, we have to look up the unique IDs for each atom
1258	>	unique_id_1 = tagRow(atom1)
1259		#else
1260	<	!! in the normal loop, the atom numbers are unique
1261	<	unique_id_1 = atom1
1260	>	!! in the normal loop, the atom numbers are unique
1261	>	unique_id_1 = atom1
1262		#endif
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	<
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	<	unique_id_2 = tagColumn(atom2)
1277	>	unique_id_2 = tagColumn(atom2)
1278		#else
1279	<	unique_id_2 = atom2
1279	>	unique_id_2 = atom2
1280		#endif
1281	<
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
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	<
1104	<	return
1105	<	end function skipThisPair
1106	<
1107	<	function FF_UsesDirectionalAtoms() result(doesit)
1108	<	logical :: doesit
1109	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1110	<	FF_uses_GB .or. FF_uses_RF
1111	<	end function FF_UsesDirectionalAtoms
1112	<
1113	<	function FF_RequiresPrepairCalc() result(doesit)
1114	<	logical :: doesit
1115	<	doesit = FF_uses_EAM
1116	<	end function FF_RequiresPrepairCalc
1117	<
1118	<	function FF_RequiresPostpairCalc() result(doesit)
1119	<	logical :: doesit
1120	<	doesit = FF_uses_RF
1121	<	end function FF_RequiresPostpairCalc
1122	<
1123	<	!! This cleans componets of force arrays belonging only to fortran
1124	<
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