[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/do

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs.
Revision 900 by chuckv, Tue Jan 6 18:54:57 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.2 2003-03-21 22:11:50 chuckv Exp $, $Date: 2003-03-21 22:11:50 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
7	>	!! @version $Id: do_Forces.F90,v 1.45 2004-01-06 18:54:57 chuckv Exp $, $Date: 2004-01-06 18:54:57 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $
8
9		module do_Forces
10		use force_globals
#	Line 17 \| Line 17 \| module do_Forces
17		use dipole_dipole
18		use reaction_field
19		use gb_pair
20	+	use vector_class
21	+	use eam
22	+	use status
23		#ifdef IS_MPI
24		use mpiSimulation
25		#endif
#	Line 27 \| Line 30 \| module do_Forces
30		#define __FORTRAN90
31		#include "fForceField.h"
32
33	<	logical, save :: do_forces_initialized = .false.
33	>	logical, save :: haveRlist = .false.
34	>	logical, save :: haveNeighborList = .false.
35	>	logical, save :: havePolicies = .false.
36	>	logical, save :: haveSIMvariables = .false.
37	>	logical, save :: havePropertyMap = .false.
38	>	logical, save :: haveSaneForceField = .false.
39		logical, save :: FF_uses_LJ
40		logical, save :: FF_uses_sticky
41		logical, save :: FF_uses_dipoles
42		logical, save :: FF_uses_RF
43		logical, save :: FF_uses_GB
44		logical, save :: FF_uses_EAM
45	+	logical, save :: SIM_uses_LJ
46	+	logical, save :: SIM_uses_sticky
47	+	logical, save :: SIM_uses_dipoles
48	+	logical, save :: SIM_uses_RF
49	+	logical, save :: SIM_uses_GB
50	+	logical, save :: SIM_uses_EAM
51	+	logical, save :: SIM_requires_postpair_calc
52	+	logical, save :: SIM_requires_prepair_calc
53	+	logical, save :: SIM_uses_directional_atoms
54	+	logical, save :: SIM_uses_PBC
55
56	+	real(kind=dp), save :: rlist, rlistsq
57	+
58		public :: init_FF
59		public :: do_force_loop
60	+	public :: setRlistDF
61
62	+
63	+	#ifdef PROFILE
64	+	public :: getforcetime
65	+	real, save :: forceTime = 0
66	+	real :: forceTimeInitial, forceTimeFinal
67	+	integer :: nLoops
68	+	#endif
69	+
70	+	type :: Properties
71	+	logical :: is_lj = .false.
72	+	logical :: is_sticky = .false.
73	+	logical :: is_dp = .false.
74	+	logical :: is_gb = .false.
75	+	logical :: is_eam = .false.
76	+	real(kind=DP) :: dipole_moment = 0.0_DP
77	+	end type Properties
78	+
79	+	type(Properties), dimension(:),allocatable :: PropertyMap
80	+
81		contains
82
83	+	subroutine setRlistDF( this_rlist )
84	+
85	+	real(kind=dp) :: this_rlist
86	+
87	+	rlist = this_rlist
88	+	rlistsq = rlist * rlist
89	+
90	+	haveRlist = .true.
91	+
92	+	end subroutine setRlistDF
93	+
94	+	subroutine createPropertyMap(status)
95	+	integer :: nAtypes
96	+	integer :: status
97	+	integer :: i
98	+	logical :: thisProperty
99	+	real (kind=DP) :: thisDPproperty
100	+
101	+	status = 0
102	+
103	+	nAtypes = getSize(atypes)
104	+
105	+	if (nAtypes == 0) then
106	+	status = -1
107	+	return
108	+	end if
109	+
110	+	if (.not. allocated(PropertyMap)) then
111	+	allocate(PropertyMap(nAtypes))
112	+	endif
113	+
114	+	do i = 1, nAtypes
115	+	call getElementProperty(atypes, i, "is_LJ", thisProperty)
116	+	PropertyMap(i)%is_LJ = thisProperty
117	+	call getElementProperty(atypes, i, "is_DP", thisProperty)
118	+	PropertyMap(i)%is_DP = thisProperty
119	+
120	+	if (thisProperty) then
121	+	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
122	+	PropertyMap(i)%dipole_moment = thisDPproperty
123	+	endif
124	+
125	+	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
126	+	PropertyMap(i)%is_Sticky = thisProperty
127	+	call getElementProperty(atypes, i, "is_GB", thisProperty)
128	+	PropertyMap(i)%is_GB = thisProperty
129	+	call getElementProperty(atypes, i, "is_EAM", thisProperty)
130	+	PropertyMap(i)%is_EAM = thisProperty
131	+	end do
132	+
133	+	havePropertyMap = .true.
134	+
135	+	end subroutine createPropertyMap
136	+
137	+	subroutine setSimVariables()
138	+	SIM_uses_LJ = SimUsesLJ()
139	+	SIM_uses_sticky = SimUsesSticky()
140	+	SIM_uses_dipoles = SimUsesDipoles()
141	+	SIM_uses_RF = SimUsesRF()
142	+	SIM_uses_GB = SimUsesGB()
143	+	SIM_uses_EAM = SimUsesEAM()
144	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
145	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
146	+	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
147	+	SIM_uses_PBC = SimUsesPBC()
148	+
149	+	haveSIMvariables = .true.
150	+
151	+	return
152	+	end subroutine setSimVariables
153	+
154	+	subroutine doReadyCheck(error)
155	+	integer, intent(out) :: error
156	+
157	+	integer :: myStatus
158	+
159	+	error = 0
160	+
161	+	if (.not. havePropertyMap) then
162	+
163	+	myStatus = 0
164	+
165	+	call createPropertyMap(myStatus)
166	+
167	+	if (myStatus .ne. 0) then
168	+	write(default_error, *) 'createPropertyMap failed in do_Forces!'
169	+	error = -1
170	+	return
171	+	endif
172	+	endif
173	+
174	+	if (.not. haveSIMvariables) then
175	+	call setSimVariables()
176	+	endif
177	+
178	+	if (.not. haveRlist) then
179	+	write(default_error, *) 'rList has not been set in do_Forces!'
180	+	error = -1
181	+	return
182	+	endif
183	+
184	+	if (SIM_uses_LJ .and. FF_uses_LJ) then
185	+	if (.not. havePolicies) then
186	+	write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
187	+	error = -1
188	+	return
189	+	endif
190	+	endif
191	+
192	+	if (.not. haveNeighborList) then
193	+	write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
194	+	error = -1
195	+	return
196	+	end if
197	+
198	+	if (.not. haveSaneForceField) then
199	+	write(default_error, *) 'Force Field is not sane in do_Forces!'
200	+	error = -1
201	+	return
202	+	end if
203	+
204	+	#ifdef IS_MPI
205	+	if (.not. isMPISimSet()) then
206	+	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
207	+	error = -1
208	+	return
209	+	endif
210	+	#endif
211	+	return
212	+	end subroutine doReadyCheck
213	+
214	+
215		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
216
217		integer, intent(in) :: LJMIXPOLICY
#	Line 84 \| Line 256 \| contains
256		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
257		if (nMatches .gt. 0) FF_uses_EAM = .true.
258
259	+	!! Assume sanity (for the sake of argument)
260	+	haveSaneForceField = .true.
261	+
262		!! check to make sure the FF_uses_RF setting makes sense
263
264		if (FF_uses_dipoles) then
90	–	rrf = getRrf()
91	–	rt = getRt()
92	–	call initialize_dipole(rrf, rt)
265		if (FF_uses_RF) then
266		dielect = getDielect()
267	<	call initialize_rf(rrf, rt, dielect)
267	>	call initialize_rf(dielect)
268		endif
269		else
270		if (FF_uses_RF) then
271		write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
272		thisStat = -1
273	+	haveSaneForceField = .false.
274		return
275		endif
276	<	endif
276	>	endif
277
278		if (FF_uses_LJ) then
279
107	–	call getRcut(rcut)
108	–
280		select case (LJMIXPOLICY)
281		case (LB_MIXING_RULE)
282	<	call init_lj_FF(LB_MIXING_RULE, rcut, my_status)
282	>	call init_lj_FF(LB_MIXING_RULE, my_status)
283		case (EXPLICIT_MIXING_RULE)
284	<	call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
284	>	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
285		case default
286		write(default_error,*) 'unknown LJ Mixing Policy!'
287		thisStat = -1
288	+	haveSaneForceField = .false.
289		return
290		end select
291		if (my_status /= 0) then
292		thisStat = -1
293	+	haveSaneForceField = .false.
294		return
295		end if
296	+	havePolicies = .true.
297		endif
298
299		if (FF_uses_sticky) then
300		call check_sticky_FF(my_status)
301	+	if (my_status /= 0) then
302	+	thisStat = -1
303	+	haveSaneForceField = .false.
304	+	return
305	+	end if
306	+	endif
307	+
308	+
309	+	if (FF_uses_EAM) then
310	+	call init_EAM_FF(my_status)
311		if (my_status /= 0) then
312	+	write(default_error, *) "init_EAM_FF returned a bad status"
313		thisStat = -1
314	+	haveSaneForceField = .false.
315		return
316		end if
317		endif
318	<
318	>
319		if (FF_uses_GB) then
320		call check_gb_pair_FF(my_status)
321		if (my_status .ne. 0) then
322		thisStat = -1
323	+	haveSaneForceField = .false.
324		return
325		endif
326		endif
327
328		if (FF_uses_GB .and. FF_uses_LJ) then
329		endif
330	<
331	<
332	<	do_forces_initialized = .true.
333	<
330	>	if (.not. haveNeighborList) then
331	>	!! Create neighbor lists
332	>	call expandNeighborList(nLocal, my_status)
333	>	if (my_Status /= 0) then
334	>	write(default_error,*) "SimSetup: ExpandNeighborList returned error."
335	>	thisStat = -1
336	>	return
337	>	endif
338	>	haveNeighborList = .true.
339	>	endif
340	>
341		end subroutine init_FF
342
343
#	Line 152 \| Line 346 \| contains
346		subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
347		error)
348		!! Position array provided by C, dimensioned by getNlocal
349	<	real ( kind = dp ), dimension(3,getNlocal()) :: q
349	>	real ( kind = dp ), dimension(3,nLocal) :: q
350		!! Rotation Matrix for each long range particle in simulation.
351	<	real( kind = dp), dimension(9,getNlocal()) :: A
351	>	real( kind = dp), dimension(9,nLocal) :: A
352		!! Unit vectors for dipoles (lab frame)
353	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
353	>	real( kind = dp ), dimension(3,nLocal) :: u_l
354		!! Force array provided by C, dimensioned by getNlocal
355	<	real ( kind = dp ), dimension(3,getNlocal()) :: f
355	>	real ( kind = dp ), dimension(3,nLocal) :: f
356		!! Torsion array provided by C, dimensioned by getNlocal
357	<	real( kind = dp ), dimension(3,getNlocal()) :: t
357	>	real( kind = dp ), dimension(3,nLocal) :: t
358	>
359		!! Stress Tensor
360		real( kind = dp), dimension(9) :: tau
361		real ( kind = dp ) :: pot
362		logical ( kind = 2) :: do_pot_c, do_stress_c
363		logical :: do_pot
364		logical :: do_stress
365	<	#ifdef IS_MPI
365	>	#ifdef IS_MPI
366		real( kind = DP ) :: pot_local
367		integer :: nrow
368		integer :: ncol
369	+	integer :: nprocs
370		#endif
175	–	integer :: nlocal
371		integer :: natoms
372		logical :: update_nlist
373		integer :: i, j, jbeg, jend, jnab
374		integer :: nlist
375	<	real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut
375	>	real( kind = DP ) :: rijsq
376		real(kind=dp),dimension(3) :: d
377		real(kind=dp) :: rfpot, mu_i, virial
378	<	integer :: me_i
378	>	integer :: me_i, me_j
379		logical :: is_dp_i
380		integer :: neighborListSize
381		integer :: listerror, error
382		integer :: localError
383	+	integer :: propPack_i, propPack_j
384
385	+	real(kind=dp) :: listSkin = 1.0
386	+
387		!! initialize local variables
388
389		#ifdef IS_MPI
390	<	nlocal = getNlocal()
390	>	pot_local = 0.0_dp
391		nrow = getNrow(plan_row)
392		ncol = getNcol(plan_col)
393		#else
196	–	nlocal = getNlocal()
394		natoms = nlocal
395		#endif
396
397	<	call getRcut(rcut,rc2=rcutsq)
201	<	call getRlist(rlist,rlistsq)
202	<
203	<	call check_initialization(localError)
397	>	call doReadyCheck(localError)
398		if ( localError .ne. 0 ) then
399	+	call handleError("do_force_loop", "Not Initialized")
400		error = -1
401		return
402		end if
#	Line 217 \| Line 412 \| contains
412		call gather(q,q_Row,plan_row3d)
413		call gather(q,q_Col,plan_col3d)
414
415	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
415	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
416		call gather(u_l,u_l_Row,plan_row3d)
417		call gather(u_l,u_l_Col,plan_col3d)
418
#	Line 226 \| Line 421 \| contains
421		endif
422
423		#endif
424	<
425	<	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
424	>
425	>	!! Begin force loop timing:
426	>	#ifdef PROFILE
427	>	call cpu_time(forceTimeInitial)
428	>	nloops = nloops + 1
429	>	#endif
430	>
431	>	if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
432		!! See if we need to update neighbor lists
433	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
433	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
434		!! if_mpi_gather_stuff_for_prepair
435		!! do_prepair_loop_if_needed
436		!! if_mpi_scatter_stuff_from_prepair
437		!! if_mpi_gather_stuff_from_prepair_to_main_loop
438	<	else
439	<	!! See if we need to update neighbor lists
440	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
438	>
439	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
440	>	#ifdef IS_MPI
441	>
442	>	if (update_nlist) then
443	>
444	>	!! save current configuration, construct neighbor list,
445	>	!! and calculate forces
446	>	call saveNeighborList(nlocal, q)
447	>
448	>	neighborListSize = size(list)
449	>	nlist = 0
450	>
451	>	do i = 1, nrow
452	>	point(i) = nlist + 1
453	>
454	>	prepair_inner: do j = 1, ncol
455	>
456	>	if (skipThisPair(i,j)) cycle prepair_inner
457	>
458	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
459	>
460	>	if (rijsq < rlistsq) then
461	>
462	>	nlist = nlist + 1
463	>
464	>	if (nlist > neighborListSize) then
465	>	call expandNeighborList(nlocal, listerror)
466	>	if (listerror /= 0) then
467	>	error = -1
468	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
469	>	return
470	>	end if
471	>	neighborListSize = size(list)
472	>	endif
473	>
474	>	list(nlist) = j
475	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
476	>	endif
477	>	enddo prepair_inner
478	>	enddo
479	>
480	>	point(nrow + 1) = nlist + 1
481	>
482	>	else !! (of update_check)
483	>
484	>	! use the list to find the neighbors
485	>	do i = 1, nrow
486	>	JBEG = POINT(i)
487	>	JEND = POINT(i+1) - 1
488	>	! check thiat molecule i has neighbors
489	>	if (jbeg .le. jend) then
490	>
491	>	do jnab = jbeg, jend
492	>	j = list(jnab)
493	>
494	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
495	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
496	>	u_l, A, f, t, pot_local)
497	>
498	>	enddo
499	>	endif
500	>	enddo
501		endif
502
503	<	#ifdef IS_MPI
503	>	#else
504
505		if (update_nlist) then
506	+
507	+	! save current configuration, contruct neighbor list,
508	+	! and calculate forces
509	+	call saveNeighborList(natoms, q)
510	+
511	+	neighborListSize = size(list)
512	+
513	+	nlist = 0
514	+
515	+	do i = 1, natoms-1
516	+	point(i) = nlist + 1
517	+
518	+	prepair_inner: do j = i+1, natoms
519	+
520	+	if (skipThisPair(i,j)) cycle prepair_inner
521	+
522	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
523	+
524	+
525	+	if (rijsq < rlistsq) then
526	+
527	+
528	+	nlist = nlist + 1
529	+
530	+	if (nlist > neighborListSize) then
531	+	call expandNeighborList(natoms, listerror)
532	+	if (listerror /= 0) then
533	+	error = -1
534	+	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
535	+	return
536	+	end if
537	+	neighborListSize = size(list)
538	+	endif
539	+
540	+	list(nlist) = j
541	+
542	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
543	+	u_l, A, f, t, pot)
544	+
545	+	endif
546	+	enddo prepair_inner
547	+	enddo
548	+
549	+	point(natoms) = nlist + 1
550
551	+	else !! (update)
552	+
553	+	! use the list to find the neighbors
554	+	do i = 1, natoms-1
555	+	JBEG = POINT(i)
556	+	JEND = POINT(i+1) - 1
557	+	! check thiat molecule i has neighbors
558	+	if (jbeg .le. jend) then
559	+
560	+	do jnab = jbeg, jend
561	+	j = list(jnab)
562	+
563	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
564	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
565	+	u_l, A, f, t, pot)
566	+
567	+	enddo
568	+	endif
569	+	enddo
570	+	endif
571	+	#endif
572	+	!! Do rest of preforce calculations
573	+	!! do necessary preforce calculations
574	+	call do_preforce(nlocal,pot)
575	+	! we have already updated the neighbor list set it to false...
576	+	update_nlist = .false.
577	+	else
578	+	!! See if we need to update neighbor lists for non pre-pair
579	+	call checkNeighborList(nlocal, q, listSkin, update_nlist)
580	+	endif
581	+
582	+
583	+
584	+
585	+
586	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
587	+
588	+
589	+
590	+
591	+
592	+	#ifdef IS_MPI
593	+
594	+	if (update_nlist) then
595		!! save current configuration, construct neighbor list,
596		!! and calculate forces
597	<	call saveNeighborList(q)
597	>	call saveNeighborList(nlocal, q)
598
599		neighborListSize = size(list)
600		nlist = 0
601
602		do i = 1, nrow
603	+
604		point(i) = nlist + 1
605
606		inner: do j = 1, ncol
#	Line 259 \| Line 609 \| contains
609
610		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
611
612	<	if (rijsq < rlistsq) then
612	>	if (rijsq < rlistsq) then
613
614		nlist = nlist + 1
615
#	Line 275 \| Line 625 \| contains
625
626		list(nlist) = j
627
628	<	if (rijsq < rcutsq) then
629	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
630	<	u_l, A, f, t,pot)
281	<	endif
628	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
629	>	u_l, A, f, t, pot_local)
630	>
631		endif
632		enddo inner
633		enddo
#	Line 299 \| Line 648 \| contains
648
649		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
650		call do_pair(i, j, rijsq, d, do_pot, do_stress, &
651	<	u_l, A, f, t,pot)
651	>	u_l, A, f, t, pot_local)
652
653		enddo
654		endif
#	Line 309 \| Line 658 \| contains
658		#else
659
660		if (update_nlist) then
661	<
661	>
662		! save current configuration, contruct neighbor list,
663		! and calculate forces
664	<	call saveNeighborList(q)
664	>	call saveNeighborList(natoms, q)
665
666		neighborListSize = size(list)
667
#	Line 328 \| Line 677 \| contains
677		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
678
679
680	<	if (rijsq < rlistsq) then
680	>	if (rijsq < rlistsq) then
681
682		nlist = nlist + 1
683
#	Line 344 \| Line 693 \| contains
693
694		list(nlist) = j
695
696	<	if (rijsq < rcutsq) then
697	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
698	<	u_l, A, f, t,pot)
350	<	endif
696	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
697	>	u_l, A, f, t, pot)
698	>
699		endif
700		enddo inner
701		enddo
#	Line 368 \| Line 716 \| contains
716
717		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
718		call do_pair(i, j, rijsq, d, do_pot, do_stress, &
719	<	u_l, A, f, t,pot)
719	>	u_l, A, f, t, pot)
720
721		enddo
722		endif
#	Line 378 \| Line 726 \| contains
726		#endif
727
728		! phew, done with main loop.
729	<
729	>
730	>	!! Do timing
731	>	#ifdef PROFILE
732	>	call cpu_time(forceTimeFinal)
733	>	forceTime = forceTime + forceTimeFinal - forceTimeInitial
734	>	#endif
735	>
736	>
737		#ifdef IS_MPI
738		!!distribute forces
739	<
740	<	call scatter(f_Row,f,plan_row3d)
739	>
740	>	f_temp = 0.0_dp
741	>	call scatter(f_Row,f_temp,plan_row3d)
742	>	do i = 1,nlocal
743	>	f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
744	>	end do
745	>
746	>	f_temp = 0.0_dp
747		call scatter(f_Col,f_temp,plan_col3d)
748		do i = 1,nlocal
749		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
750		end do
751
752	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
753	<	call scatter(t_Row,t,plan_row3d)
752	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
753	>	t_temp = 0.0_dp
754	>	call scatter(t_Row,t_temp,plan_row3d)
755	>	do i = 1,nlocal
756	>	t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
757	>	end do
758	>	t_temp = 0.0_dp
759		call scatter(t_Col,t_temp,plan_col3d)
760
761		do i = 1,nlocal
#	Line 400 \| Line 766 \| contains
766		if (do_pot) then
767		! scatter/gather pot_row into the members of my column
768		call scatter(pot_Row, pot_Temp, plan_row)
769	<
769	>
770		! scatter/gather pot_local into all other procs
771		! add resultant to get total pot
772		do i = 1, nlocal
773		pot_local = pot_local + pot_Temp(i)
774		enddo
775	+
776	+	pot_Temp = 0.0_DP
777
410	–	pot_Temp = 0.0_DP
411	–
778		call scatter(pot_Col, pot_Temp, plan_col)
779		do i = 1, nlocal
780		pot_local = pot_local + pot_Temp(i)
781		enddo
782	<
782	>
783		endif
784		#endif
785
786	<	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
786	>	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
787
788	<	if (FF_uses_RF .and. SimUsesRF()) then
788	>	if (FF_uses_RF .and. SIM_uses_RF) then
789
790		#ifdef IS_MPI
791		call scatter(rf_Row,rf,plan_row3d)
#	Line 429 \| Line 795 \| contains
795		end do
796		#endif
797
798	<	do i = 1, getNlocal()
798	>	do i = 1, nLocal
799
800		rfpot = 0.0_DP
801		#ifdef IS_MPI
#	Line 437 \| Line 803 \| contains
803		#else
804		me_i = atid(i)
805		#endif
806	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
807	<	if ( is_DP_i ) then
808	<	call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
806	>
807	>	if (PropertyMap(me_i)%is_DP) then
808	>
809	>	mu_i = PropertyMap(me_i)%dipole_moment
810	>
811		!! The reaction field needs to include a self contribution
812		!! to the field:
813	<	call accumulate_self_rf(i, mu_i, u_l)
813	>	call accumulate_self_rf(i, mu_i, u_l)
814		!! Get the reaction field contribution to the
815		!! potential and torques:
816		call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#	Line 461 \| Line 829 \| contains
829		#ifdef IS_MPI
830
831		if (do_pot) then
832	<	pot = pot_local
832	>	pot = pot + pot_local
833		!! we assume the c code will do the allreduce to get the total potential
834		!! we could do it right here if we needed to...
835		endif
836
837		if (do_stress) then
838	<	call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
838	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
839		mpi_comm_world,mpi_err)
840	<	call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
840	>	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
841		mpi_comm_world,mpi_err)
842		endif
843
#	Line 479 \| Line 847 \| contains
847		tau = tau_Temp
848		virial = virial_Temp
849		endif
850	<
850	>
851		#endif
852
853	+
854	+
855		end subroutine do_force_loop
856
857	<	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t,pot)
857	>	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
858
859		real( kind = dp ) :: pot
860	<	real( kind = dp ), dimension(:,:) :: u_l
861	<	real (kind=dp), dimension(:,:) :: A
862	<	real (kind=dp), dimension(:,:) :: f
863	<	real (kind=dp), dimension(:,:) :: t
860	>	real( kind = dp ), dimension(3,nLocal) :: u_l
861	>	real (kind=dp), dimension(9,nLocal) :: A
862	>	real (kind=dp), dimension(3,nLocal) :: f
863	>	real (kind=dp), dimension(3,nLocal) :: t
864
865		logical, intent(inout) :: do_pot, do_stress
866		integer, intent(in) :: i, j
#	Line 500 \| Line 870 \| contains
870		logical :: is_LJ_i, is_LJ_j
871		logical :: is_DP_i, is_DP_j
872		logical :: is_GB_i, is_GB_j
873	+	logical :: is_EAM_i,is_EAM_j
874		logical :: is_Sticky_i, is_Sticky_j
875		integer :: me_i, me_j
876	<
876	>	integer :: propPack_i
877	>	integer :: propPack_j
878		r = sqrt(rijsq)
879
880		#ifdef IS_MPI
881	+	if (tagRow(i) .eq. tagColumn(j)) then
882	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
883	+	endif
884
885		me_i = atid_row(i)
886		me_j = atid_col(j)
#	Line 516 \| Line 891 \| contains
891		me_j = atid(j)
892
893		#endif
894	+
895	+	if (FF_uses_LJ .and. SIM_uses_LJ) then
896
897	<	if (FF_uses_LJ .and. SimUsesLJ()) then
521	<	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
522	<	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
523	<
524	<	if ( is_LJ_i .and. is_LJ_j ) &
897	>	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) &
898		call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
899	+
900		endif
901
902	<	if (FF_uses_dipoles .and. SimUsesDipoles()) then
529	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
530	<	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
902	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
903
904	<	if ( is_DP_i .and. is_DP_j ) then
533	<
904	>	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
905		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
906		do_pot, do_stress)
907	<	if (FF_uses_RF .and. SimUsesRF()) then
907	>	if (FF_uses_RF .and. SIM_uses_RF) then
908		call accumulate_rf(i, j, r, u_l)
909		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
910		endif
#	Line 541 \| Line 912 \| contains
912		endif
913		endif
914
915	<	if (FF_uses_Sticky .and. SimUsesSticky()) then
915	>	if (FF_uses_Sticky .and. SIM_uses_sticky) then
916
917	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
547	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
548	<
549	<	if ( is_Sticky_i .and. is_Sticky_j ) then
917	>	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
918		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
919		do_pot, do_stress)
920		endif
921		endif
922
923
924	<	if (FF_uses_GB .and. SimUsesGB()) then
557	<
558	<	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
559	<	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
924	>	if (FF_uses_GB .and. SIM_uses_GB) then
925
926	<	if ( is_GB_i .and. is_GB_j ) then
926	>	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
927		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
928		do_pot, do_stress)
929		endif
930	+
931		endif
932
933	+
934	+
935	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
936	+
937	+	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
938	+	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
939	+	endif
940	+
941	+	endif
942	+
943		end subroutine do_pair
944
945
946	+
947	+	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
948	+	real( kind = dp ) :: pot
949	+	real( kind = dp ), dimension(3,nLocal) :: u_l
950	+	real (kind=dp), dimension(9,nLocal) :: A
951	+	real (kind=dp), dimension(3,nLocal) :: f
952	+	real (kind=dp), dimension(3,nLocal) :: t
953	+
954	+	logical, intent(inout) :: do_pot, do_stress
955	+	integer, intent(in) :: i, j
956	+	real ( kind = dp ), intent(inout) :: rijsq
957	+	real ( kind = dp ) :: r
958	+	real ( kind = dp ), intent(inout) :: d(3)
959	+
960	+	logical :: is_EAM_i, is_EAM_j
961	+
962	+	integer :: me_i, me_j
963	+
964	+	r = sqrt(rijsq)
965	+
966	+
967	+	#ifdef IS_MPI
968	+	if (tagRow(i) .eq. tagColumn(j)) then
969	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
970	+	endif
971	+
972	+	me_i = atid_row(i)
973	+	me_j = atid_col(j)
974	+
975	+	#else
976	+
977	+	me_i = atid(i)
978	+	me_j = atid(j)
979	+
980	+	#endif
981	+
982	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
983	+
984	+	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
985	+	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
986	+
987	+	endif
988	+
989	+	end subroutine do_prepair
990	+
991	+
992	+
993	+
994	+	subroutine do_preforce(nlocal,pot)
995	+	integer :: nlocal
996	+	real( kind = dp ) :: pot
997	+
998	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
999	+	call calc_EAM_preforce_Frho(nlocal,pot)
1000	+	endif
1001	+
1002	+
1003	+	end subroutine do_preforce
1004	+
1005	+
1006		subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1007
1008		real (kind = dp), dimension(3) :: q_i
1009		real (kind = dp), dimension(3) :: q_j
1010		real ( kind = dp ), intent(out) :: r_sq
1011	<	real( kind = dp ) :: d(3)
1012	<	real( kind = dp ) :: d_old(3)
1013	<	d(1:3) = q_i(1:3) - q_j(1:3)
1014	<	d_old = d
1011	>	real( kind = dp ) :: d(3), scaled(3)
1012	>	integer i
1013	>
1014	>	d(1:3) = q_j(1:3) - q_i(1:3)
1015	>
1016		! Wrap back into periodic box if necessary
1017	<	if ( SimUsesPBC() ) then
1017	>	if ( SIM_uses_PBC ) then
1018	>
1019	>	if( .not.boxIsOrthorhombic ) then
1020	>	! calc the scaled coordinates.
1021	>
1022	>	scaled = matmul(HmatInv, d)
1023	>
1024	>	! wrap the scaled coordinates
1025
1026	<	d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
1027	<	int(abs(d(1:3)/box(1:3)) + 0.5_dp)
1026	>	scaled = scaled - anint(scaled)
1027	>
1028
1029	<	endif
1030	<	r_sq = dot_product(d,d)
587	<
588	<	end subroutine get_interatomic_vector
1029	>	! calc the wrapped real coordinates from the wrapped scaled
1030	>	! coordinates
1031
1032	<	subroutine check_initialization(error)
591	<	integer, intent(out) :: error
592	<
593	<	error = 0
594	<	! Make sure we are properly initialized.
595	<	if (.not. do_forces_initialized) then
596	<	error = -1
597	<	return
598	<	endif
1032	>	d = matmul(Hmat,scaled)
1033
1034	<	#ifdef IS_MPI
1035	<	if (.not. isMPISimSet()) then
1036	<	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
1037	<	error = -1
1038	<	return
1034	>	else
1035	>	! calc the scaled coordinates.
1036	>
1037	>	do i = 1, 3
1038	>	scaled(i) = d(i) * HmatInv(i,i)
1039	>
1040	>	! wrap the scaled coordinates
1041	>
1042	>	scaled(i) = scaled(i) - anint(scaled(i))
1043	>
1044	>	! calc the wrapped real coordinates from the wrapped scaled
1045	>	! coordinates
1046	>
1047	>	d(i) = scaled(i)*Hmat(i,i)
1048	>	enddo
1049	>	endif
1050	>
1051		endif
606	–	#endif
1052
1053	<	return
1054	<	end subroutine check_initialization
1055	<
1053	>	r_sq = dot_product(d,d)
1054	>
1055	>	end subroutine get_interatomic_vector
1056
1057		subroutine zero_work_arrays()
1058
#	Line 640 \| Line 1085 \| contains
1085
1086		#endif
1087
1088	+
1089	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
1090	+	call clean_EAM()
1091	+	endif
1092	+
1093	+
1094	+
1095	+
1096	+
1097		rf = 0.0_dp
1098		tau_Temp = 0.0_dp
1099		virial_Temp = 0.0_dp
646	–
1100		end subroutine zero_work_arrays
1101
1102		function skipThisPair(atom1, atom2) result(skip_it)
#	Line 687 \| Line 1140 \| contains
1140		#else
1141		unique_id_2 = atom2
1142		#endif
1143	<
1143	>
1144		#ifdef IS_MPI
1145		!! this situation should only arise in MPI simulations
1146		if (unique_id_1 == unique_id_2) then
#	Line 697 \| Line 1150 \| contains
1150
1151		!! this prevents us from doing the pair on multiple processors
1152		if (unique_id_1 < unique_id_2) then
1153	<	if (mod(unique_id_1 + unique_id_2,2) == 0) skip_it = .true.
1154	<	return
1153	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1154	>	skip_it = .true.
1155	>	return
1156	>	endif
1157		else
1158	<	if (mod(unique_id_1 + unique_id_2,2) == 1) skip_it = .true.
1159	<	return
1158	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1159	>	skip_it = .true.
1160	>	return
1161	>	endif
1162		endif
1163		#endif
1164	<
1164	>
1165		!! the rest of these situations can happen in all simulations:
1166		do i = 1, nExcludes_global
1167		if ((excludesGlobal(i) == unique_id_1) .or. &
#	Line 713 \| Line 1170 \| contains
1170		return
1171		endif
1172		enddo
1173	<
1173	>
1174		do i = 1, nExcludes_local
1175		if (excludesLocal(1,i) == unique_id_1) then
1176		if (excludesLocal(2,i) == unique_id_2) then
#	Line 749 \| Line 1206 \| end module do_Forces
1206		doesit = FF_uses_RF
1207		end function FF_RequiresPostpairCalc
1208
1209	+	#ifdef PROFILE
1210	+	function getforcetime() result(totalforcetime)
1211	+	real(kind=dp) :: totalforcetime
1212	+	totalforcetime = forcetime
1213	+	end function getforcetime
1214	+	#endif
1215	+
1216	+	!! This cleans componets of force arrays belonging only to fortran
1217	+
1218		end module do_Forces

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents): Revision 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs. Revision 900 by chuckv, Tue Jan 6 18:54:57 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs.
Revision 900 by chuckv, Tue Jan 6 18:54:57 2004 UTC