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 872 by chrisfen, Fri Nov 21 19:31:05 2003 UTC vs.
Revision 1144 by tim, Sat May 1 18:52:38 2004 UTC

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

Diff Legend

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