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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC vs.
Revision 1132 by tim, Sat Apr 24 04:31:36 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.15 2003-04-11 15:16:59 gezelter Exp $, $Date: 2003-04-11 15:16:59 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $
7	>	!! @version $Id: do_Forces.F90,v 1.51 2004-04-24 04:31:36 tim Exp $, $Date: 2004-04-24 04:31:36 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
8
9		module do_Forces
10		use force_globals
#	Line 15 \| Line 15 \| module do_Forces
15		use lj
16		use sticky_pair
17		use dipole_dipole
18	+	use charge_charge
19		use reaction_field
20		use gb_pair
21	+	use vector_class
22	+	use eam
23	+	use status
24		#ifdef IS_MPI
25		use mpiSimulation
26		#endif
#	Line 27 \| Line 31 \| module do_Forces
31		#define __FORTRAN90
32		#include "fForceField.h"
33
34	<	logical, save :: do_forces_initialized = .false.
34	>	logical, save :: haveRlist = .false.
35	>	logical, save :: haveNeighborList = .false.
36	>	logical, save :: havePolicies = .false.
37	>	logical, save :: haveSIMvariables = .false.
38	>	logical, save :: havePropertyMap = .false.
39	>	logical, save :: haveSaneForceField = .false.
40		logical, save :: FF_uses_LJ
41		logical, save :: FF_uses_sticky
42	+	logical, save :: FF_uses_charges
43		logical, save :: FF_uses_dipoles
44		logical, save :: FF_uses_RF
45		logical, save :: FF_uses_GB
46		logical, save :: FF_uses_EAM
47	+	logical, save :: SIM_uses_LJ
48	+	logical, save :: SIM_uses_sticky
49	+	logical, save :: SIM_uses_charges
50	+	logical, save :: SIM_uses_dipoles
51	+	logical, save :: SIM_uses_RF
52	+	logical, save :: SIM_uses_GB
53	+	logical, save :: SIM_uses_EAM
54	+	logical, save :: SIM_requires_postpair_calc
55	+	logical, save :: SIM_requires_prepair_calc
56	+	logical, save :: SIM_uses_directional_atoms
57	+	logical, save :: SIM_uses_PBC
58
59	+	real(kind=dp), save :: rlist, rlistsq
60	+
61		public :: init_FF
62		public :: do_force_loop
63	+	public :: setRlistDF
64
65	+
66	+	#ifdef PROFILE
67	+	public :: getforcetime
68	+	real, save :: forceTime = 0
69	+	real :: forceTimeInitial, forceTimeFinal
70	+	integer :: nLoops
71	+	#endif
72	+
73	+	type :: Properties
74	+	logical :: is_lj = .false.
75	+	logical :: is_sticky = .false.
76	+	logical :: is_dp = .false.
77	+	logical :: is_gb = .false.
78	+	logical :: is_eam = .false.
79	+	logical :: is_charge = .false.
80	+	real(kind=DP) :: charge = 0.0_DP
81	+	real(kind=DP) :: dipole_moment = 0.0_DP
82	+	end type Properties
83	+
84	+	type(Properties), dimension(:),allocatable :: PropertyMap
85	+
86		contains
87
88	+	subroutine setRlistDF( this_rlist )
89	+
90	+	real(kind=dp) :: this_rlist
91	+
92	+	rlist = this_rlist
93	+	rlistsq = rlist * rlist
94	+
95	+	haveRlist = .true.
96	+
97	+	end subroutine setRlistDF
98	+
99	+	subroutine createPropertyMap(status)
100	+	integer :: nAtypes
101	+	integer :: status
102	+	integer :: i
103	+	logical :: thisProperty
104	+	real (kind=DP) :: thisDPproperty
105	+
106	+	status = 0
107	+
108	+	nAtypes = getSize(atypes)
109	+
110	+	if (nAtypes == 0) then
111	+	status = -1
112	+	return
113	+	end if
114	+
115	+	if (.not. allocated(PropertyMap)) then
116	+	allocate(PropertyMap(nAtypes))
117	+	endif
118	+
119	+	do i = 1, nAtypes
120	+	call getElementProperty(atypes, i, "is_LJ", thisProperty)
121	+	PropertyMap(i)%is_LJ = thisProperty
122	+
123	+	call getElementProperty(atypes, i, "is_Charge", thisProperty)
124	+	PropertyMap(i)%is_Charge = thisProperty
125	+
126	+	if (thisProperty) then
127	+	call getElementProperty(atypes, i, "charge", thisDPproperty)
128	+	PropertyMap(i)%charge = thisDPproperty
129	+	endif
130	+
131	+	call getElementProperty(atypes, i, "is_DP", thisProperty)
132	+	PropertyMap(i)%is_DP = thisProperty
133	+
134	+	if (thisProperty) then
135	+	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
136	+	PropertyMap(i)%dipole_moment = thisDPproperty
137	+	endif
138	+
139	+	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
140	+	PropertyMap(i)%is_Sticky = thisProperty
141	+	call getElementProperty(atypes, i, "is_GB", thisProperty)
142	+	PropertyMap(i)%is_GB = thisProperty
143	+	call getElementProperty(atypes, i, "is_EAM", thisProperty)
144	+	PropertyMap(i)%is_EAM = thisProperty
145	+	end do
146	+
147	+	havePropertyMap = .true.
148	+
149	+	end subroutine createPropertyMap
150	+
151	+	subroutine setSimVariables()
152	+	SIM_uses_LJ = SimUsesLJ()
153	+	SIM_uses_sticky = SimUsesSticky()
154	+	SIM_uses_charges = SimUsesCharges()
155	+	SIM_uses_dipoles = SimUsesDipoles()
156	+	SIM_uses_RF = SimUsesRF()
157	+	SIM_uses_GB = SimUsesGB()
158	+	SIM_uses_EAM = SimUsesEAM()
159	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
160	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
161	+	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
162	+	SIM_uses_PBC = SimUsesPBC()
163	+
164	+	haveSIMvariables = .true.
165	+
166	+	return
167	+	end subroutine setSimVariables
168	+
169	+	subroutine doReadyCheck(error)
170	+	integer, intent(out) :: error
171	+
172	+	integer :: myStatus
173	+
174	+	error = 0
175	+
176	+	if (.not. havePropertyMap) then
177	+
178	+	myStatus = 0
179	+
180	+	call createPropertyMap(myStatus)
181	+
182	+	if (myStatus .ne. 0) then
183	+	write(default_error, *) 'createPropertyMap failed in do_Forces!'
184	+	error = -1
185	+	return
186	+	endif
187	+	endif
188	+
189	+	if (.not. haveSIMvariables) then
190	+	call setSimVariables()
191	+	endif
192	+
193	+	if (.not. haveRlist) then
194	+	write(default_error, *) 'rList has not been set in do_Forces!'
195	+	error = -1
196	+	return
197	+	endif
198	+
199	+	if (SIM_uses_LJ .and. FF_uses_LJ) then
200	+	if (.not. havePolicies) then
201	+	write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
202	+	error = -1
203	+	return
204	+	endif
205	+	endif
206	+
207	+	if (.not. haveNeighborList) then
208	+	write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
209	+	error = -1
210	+	return
211	+	end if
212	+
213	+	if (.not. haveSaneForceField) then
214	+	write(default_error, *) 'Force Field is not sane in do_Forces!'
215	+	error = -1
216	+	return
217	+	end if
218	+
219	+	#ifdef IS_MPI
220	+	if (.not. isMPISimSet()) then
221	+	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
222	+	error = -1
223	+	return
224	+	endif
225	+	#endif
226	+	return
227	+	end subroutine doReadyCheck
228	+
229	+
230		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
231
232		integer, intent(in) :: LJMIXPOLICY
#	Line 64 \| Line 251 \| contains
251
252		FF_uses_LJ = .false.
253		FF_uses_sticky = .false.
254	+	FF_uses_charges = .false.
255		FF_uses_dipoles = .false.
256		FF_uses_GB = .false.
257		FF_uses_EAM = .false.
258
259		call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
260		if (nMatches .gt. 0) FF_uses_LJ = .true.
261	<
261	>
262	>	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
263	>	if (nMatches .gt. 0) FF_uses_charges = .true.
264	>
265		call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
266		if (nMatches .gt. 0) FF_uses_dipoles = .true.
267
#	Line 84 \| Line 275 \| contains
275		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
276		if (nMatches .gt. 0) FF_uses_EAM = .true.
277
278	+	!! Assume sanity (for the sake of argument)
279	+	haveSaneForceField = .true.
280	+	!!
281	+	if (FF_uses_charges) then
282	+	dielect = getDielect()
283	+	call initialize_charge(dielect)
284	+	endif
285	+
286	+
287		!! check to make sure the FF_uses_RF setting makes sense
288
289		if (FF_uses_dipoles) then
90	–	rrf = getRrf()
91	–	rt = getRt()
92	–	call initialize_dipole(rrf, rt)
290		if (FF_uses_RF) then
291		dielect = getDielect()
292	<	call initialize_rf(rrf, rt, dielect)
292	>	call initialize_rf(dielect)
293		endif
294		else
295		if (FF_uses_RF) then
296		write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
297		thisStat = -1
298	+	haveSaneForceField = .false.
299		return
300		endif
301	<	endif
301	>	endif
302
303		if (FF_uses_LJ) then
304
107	–	call getRcut(rcut)
108	–
305		select case (LJMIXPOLICY)
306		case (LB_MIXING_RULE)
307	<	call init_lj_FF(LB_MIXING_RULE, rcut, my_status)
307	>	call init_lj_FF(LB_MIXING_RULE, my_status)
308		case (EXPLICIT_MIXING_RULE)
309	<	call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
309	>	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
310		case default
311		write(default_error,*) 'unknown LJ Mixing Policy!'
312		thisStat = -1
313	+	haveSaneForceField = .false.
314		return
315		end select
316		if (my_status /= 0) then
317		thisStat = -1
318	+	haveSaneForceField = .false.
319		return
320		end if
321	+	havePolicies = .true.
322		endif
323
324		if (FF_uses_sticky) then
325		call check_sticky_FF(my_status)
326		if (my_status /= 0) then
327		thisStat = -1
328	+	haveSaneForceField = .false.
329		return
330		end if
331		endif
332	<
332	>
333	>
334	>	if (FF_uses_EAM) then
335	>	call init_EAM_FF(my_status)
336	>	if (my_status /= 0) then
337	>	write(default_error, *) "init_EAM_FF returned a bad status"
338	>	thisStat = -1
339	>	haveSaneForceField = .false.
340	>	return
341	>	end if
342	>	endif
343	>
344		if (FF_uses_GB) then
345		call check_gb_pair_FF(my_status)
346		if (my_status .ne. 0) then
347		thisStat = -1
348	+	haveSaneForceField = .false.
349		return
350		endif
351		endif
352
353		if (FF_uses_GB .and. FF_uses_LJ) then
354		endif
355	<	if (.not. do_forces_initialized) then
355	>	if (.not. haveNeighborList) then
356		!! Create neighbor lists
357	<	call expandNeighborList(getNlocal(), my_status)
357	>	call expandNeighborList(nLocal, my_status)
358		if (my_Status /= 0) then
359		write(default_error,*) "SimSetup: ExpandNeighborList returned error."
360		thisStat = -1
361		return
362		endif
363	+	haveNeighborList = .true.
364		endif
365	<
153	<	do_forces_initialized = .true.
154	<
365	>
366		end subroutine init_FF
367
368
#	Line 160 \| Line 371 \| contains
371		subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
372		error)
373		!! Position array provided by C, dimensioned by getNlocal
374	<	real ( kind = dp ), dimension(3,getNlocal()) :: q
374	>	real ( kind = dp ), dimension(3,nLocal) :: q
375		!! Rotation Matrix for each long range particle in simulation.
376	<	real( kind = dp), dimension(9,getNlocal()) :: A
376	>	real( kind = dp), dimension(9,nLocal) :: A
377		!! Unit vectors for dipoles (lab frame)
378	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
378	>	real( kind = dp ), dimension(3,nLocal) :: u_l
379		!! Force array provided by C, dimensioned by getNlocal
380	<	real ( kind = dp ), dimension(3,getNlocal()) :: f
380	>	real ( kind = dp ), dimension(3,nLocal) :: f
381		!! Torsion array provided by C, dimensioned by getNlocal
382	<	real( kind = dp ), dimension(3,getNlocal()) :: t
382	>	real( kind = dp ), dimension(3,nLocal) :: t
383	>
384		!! Stress Tensor
385		real( kind = dp), dimension(9) :: tau
386		real ( kind = dp ) :: pot
#	Line 179 \| Line 391 \| contains
391		real( kind = DP ) :: pot_local
392		integer :: nrow
393		integer :: ncol
394	+	integer :: nprocs
395		#endif
183	–	integer :: nlocal
396		integer :: natoms
397		logical :: update_nlist
398		integer :: i, j, jbeg, jend, jnab
399		integer :: nlist
400	<	real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut
400	>	real( kind = DP ) :: rijsq
401		real(kind=dp),dimension(3) :: d
402		real(kind=dp) :: rfpot, mu_i, virial
403	<	integer :: me_i
403	>	integer :: me_i, me_j
404		logical :: is_dp_i
405		integer :: neighborListSize
406		integer :: listerror, error
407		integer :: localError
408	+	integer :: propPack_i, propPack_j
409
410	+	real(kind=dp) :: listSkin = 1.0
411	+
412		!! initialize local variables
413
414		#ifdef IS_MPI
415		pot_local = 0.0_dp
201	–	nlocal = getNlocal()
416		nrow = getNrow(plan_row)
417		ncol = getNcol(plan_col)
418		#else
205	–	nlocal = getNlocal()
419		natoms = nlocal
420		#endif
421	<
422	<	call getRcut(rcut,rc2=rcutsq)
210	<	call getRlist(rlist,rlistsq)
211	<
212	<	call check_initialization(localError)
421	>
422	>	call doReadyCheck(localError)
423		if ( localError .ne. 0 ) then
424	+	call handleError("do_force_loop", "Not Initialized")
425		error = -1
426		return
427		end if
428		call zero_work_arrays()
429
430	+
431		do_pot = do_pot_c
432		do_stress = do_stress_c
433
#	Line 226 \| Line 438 \| contains
438		call gather(q,q_Row,plan_row3d)
439		call gather(q,q_Col,plan_col3d)
440
441	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
441	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
442		call gather(u_l,u_l_Row,plan_row3d)
443		call gather(u_l,u_l_Col,plan_col3d)
444
#	Line 235 \| Line 447 \| contains
447		endif
448
449		#endif
450	<
451	<	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
450	>
451	>	!! Begin force loop timing:
452	>	#ifdef PROFILE
453	>	call cpu_time(forceTimeInitial)
454	>	nloops = nloops + 1
455	>	#endif
456	>
457	>	if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
458		!! See if we need to update neighbor lists
459	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
459	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
460		!! if_mpi_gather_stuff_for_prepair
461		!! do_prepair_loop_if_needed
462		!! if_mpi_scatter_stuff_from_prepair
463		!! if_mpi_gather_stuff_from_prepair_to_main_loop
464	<	else
465	<	!! See if we need to update neighbor lists
466	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
464	>
465	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
466	>	#ifdef IS_MPI
467	>
468	>	if (update_nlist) then
469	>
470	>	!! save current configuration, construct neighbor list,
471	>	!! and calculate forces
472	>	call saveNeighborList(nlocal, q)
473	>
474	>	neighborListSize = size(list)
475	>	nlist = 0
476	>
477	>	do i = 1, nrow
478	>	point(i) = nlist + 1
479	>
480	>	prepair_inner: do j = 1, ncol
481	>
482	>	if (skipThisPair(i,j)) cycle prepair_inner
483	>
484	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
485	>
486	>	if (rijsq < rlistsq) then
487	>
488	>	nlist = nlist + 1
489	>
490	>	if (nlist > neighborListSize) then
491	>	call expandNeighborList(nlocal, listerror)
492	>	if (listerror /= 0) then
493	>	error = -1
494	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
495	>	return
496	>	end if
497	>	neighborListSize = size(list)
498	>	endif
499	>
500	>	list(nlist) = j
501	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
502	>	endif
503	>	enddo prepair_inner
504	>	enddo
505	>
506	>	point(nrow + 1) = nlist + 1
507	>
508	>	else !! (of update_check)
509	>
510	>	! use the list to find the neighbors
511	>	do i = 1, nrow
512	>	JBEG = POINT(i)
513	>	JEND = POINT(i+1) - 1
514	>	! check thiat molecule i has neighbors
515	>	if (jbeg .le. jend) then
516	>
517	>	do jnab = jbeg, jend
518	>	j = list(jnab)
519	>
520	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
521	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
522	>	u_l, A, f, t, pot_local)
523	>
524	>	enddo
525	>	endif
526	>	enddo
527		endif
528
529	<	#ifdef IS_MPI
529	>	#else
530
531		if (update_nlist) then
532	+
533	+	! save current configuration, contruct neighbor list,
534	+	! and calculate forces
535	+	call saveNeighborList(natoms, q)
536	+
537	+	neighborListSize = size(list)
538	+
539	+	nlist = 0
540	+
541	+	do i = 1, natoms-1
542	+	point(i) = nlist + 1
543	+
544	+	prepair_inner: do j = i+1, natoms
545	+
546	+	if (skipThisPair(i,j)) cycle prepair_inner
547	+
548	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
549	+
550	+
551	+	if (rijsq < rlistsq) then
552	+
553	+
554	+	nlist = nlist + 1
555	+
556	+	if (nlist > neighborListSize) then
557	+	call expandNeighborList(natoms, listerror)
558	+	if (listerror /= 0) then
559	+	error = -1
560	+	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
561	+	return
562	+	end if
563	+	neighborListSize = size(list)
564	+	endif
565	+
566	+	list(nlist) = j
567	+
568	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
569	+	u_l, A, f, t, pot)
570	+
571	+	endif
572	+	enddo prepair_inner
573	+	enddo
574
575	+	point(natoms) = nlist + 1
576	+
577	+	else !! (update)
578	+
579	+	! use the list to find the neighbors
580	+	do i = 1, natoms-1
581	+	JBEG = POINT(i)
582	+	JEND = POINT(i+1) - 1
583	+	! check thiat molecule i has neighbors
584	+	if (jbeg .le. jend) then
585	+
586	+	do jnab = jbeg, jend
587	+	j = list(jnab)
588	+
589	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
590	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
591	+	u_l, A, f, t, pot)
592	+
593	+	enddo
594	+	endif
595	+	enddo
596	+	endif
597	+	#endif
598	+	!! Do rest of preforce calculations
599	+	!! do necessary preforce calculations
600	+	call do_preforce(nlocal,pot)
601	+	! we have already updated the neighbor list set it to false...
602	+	update_nlist = .false.
603	+	else
604	+	!! See if we need to update neighbor lists for non pre-pair
605	+	call checkNeighborList(nlocal, q, listSkin, update_nlist)
606	+	endif
607	+
608	+
609	+
610	+
611	+
612	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
613	+
614	+
615	+
616	+
617	+
618	+	#ifdef IS_MPI
619	+
620	+	if (update_nlist) then
621		!! save current configuration, construct neighbor list,
622		!! and calculate forces
623		call saveNeighborList(nlocal, q)
#	Line 260 \| Line 626 \| contains
626		nlist = 0
627
628		do i = 1, nrow
629	+
630		point(i) = nlist + 1
631
632		inner: do j = 1, ncol
#	Line 268 \| Line 635 \| contains
635
636		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
637
638	<	if (rijsq < rlistsq) then
638	>	if (rijsq < rlistsq) then
639
640		nlist = nlist + 1
641
#	Line 284 \| Line 651 \| contains
651
652		list(nlist) = j
653
654	<	if (rijsq < rcutsq) then
655	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
656	<	u_l, A, f, t, pot_local)
290	<	endif
654	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
655	>	u_l, A, f, t, pot_local)
656	>
657		endif
658		enddo inner
659		enddo
#	Line 318 \| Line 684 \| contains
684		#else
685
686		if (update_nlist) then
687	<
687	>
688		! save current configuration, contruct neighbor list,
689		! and calculate forces
690		call saveNeighborList(natoms, q)
#	Line 337 \| Line 703 \| contains
703		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
704
705
706	<	if (rijsq < rlistsq) then
706	>	if (rijsq < rlistsq) then
707
708		nlist = nlist + 1
709
#	Line 353 \| Line 719 \| contains
719
720		list(nlist) = j
721
722	<	if (rijsq < rcutsq) then
357	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
722	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
723		u_l, A, f, t, pot)
724	<	endif
724	>
725		endif
726		enddo inner
727		enddo
#	Line 387 \| Line 752 \| contains
752		#endif
753
754		! phew, done with main loop.
755	<
755	>
756	>	!! Do timing
757	>	#ifdef PROFILE
758	>	call cpu_time(forceTimeFinal)
759	>	forceTime = forceTime + forceTimeFinal - forceTimeInitial
760	>	#endif
761	>
762	>
763		#ifdef IS_MPI
764		!!distribute forces
765
#	Line 403 \| Line 775 \| contains
775		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
776		end do
777
778	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
778	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
779		t_temp = 0.0_dp
780		call scatter(t_Row,t_temp,plan_row3d)
781		do i = 1,nlocal
#	Line 437 \| Line 809 \| contains
809		endif
810		#endif
811
812	<	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
812	>	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
813
814	<	if (FF_uses_RF .and. SimUsesRF()) then
815	<
814	>	if (FF_uses_RF .and. SIM_uses_RF) then
815	>
816		#ifdef IS_MPI
817		call scatter(rf_Row,rf,plan_row3d)
818		call scatter(rf_Col,rf_Temp,plan_col3d)
#	Line 449 \| Line 821 \| contains
821		end do
822		#endif
823
824	<	do i = 1, getNlocal()
824	>	do i = 1, nLocal
825
826		rfpot = 0.0_DP
827		#ifdef IS_MPI
#	Line 457 \| Line 829 \| contains
829		#else
830		me_i = atid(i)
831		#endif
832	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
833	<	if ( is_DP_i ) then
834	<	call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
832	>
833	>	if (PropertyMap(me_i)%is_DP) then
834	>
835	>	mu_i = PropertyMap(me_i)%dipole_moment
836	>
837		!! The reaction field needs to include a self contribution
838		!! to the field:
839	<	call accumulate_self_rf(i, mu_i, u_l)
839	>	call accumulate_self_rf(i, mu_i, u_l)
840		!! Get the reaction field contribution to the
841		!! potential and torques:
842		call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#	Line 499 \| Line 873 \| contains
873		tau = tau_Temp
874		virial = virial_Temp
875		endif
876	<
876	>
877		#endif
878
879	+
880		end subroutine do_force_loop
881
882		subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
883
884		real( kind = dp ) :: pot
885	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
886	<	real (kind=dp), dimension(9,getNlocal()) :: A
887	<	real (kind=dp), dimension(3,getNlocal()) :: f
888	<	real (kind=dp), dimension(3,getNlocal()) :: t
885	>	real( kind = dp ), dimension(3,nLocal) :: u_l
886	>	real (kind=dp), dimension(9,nLocal) :: A
887	>	real (kind=dp), dimension(3,nLocal) :: f
888	>	real (kind=dp), dimension(3,nLocal) :: t
889
890		logical, intent(inout) :: do_pot, do_stress
891		integer, intent(in) :: i, j
892		real ( kind = dp ), intent(inout) :: rijsq
893		real ( kind = dp ) :: r
894		real ( kind = dp ), intent(inout) :: d(3)
520	–	logical :: is_LJ_i, is_LJ_j
521	–	logical :: is_DP_i, is_DP_j
522	–	logical :: is_GB_i, is_GB_j
523	–	logical :: is_Sticky_i, is_Sticky_j
895		integer :: me_i, me_j
896
897	+
898		r = sqrt(rijsq)
899
528	–
529	–
900		#ifdef IS_MPI
901		if (tagRow(i) .eq. tagColumn(j)) then
902		write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
903		endif
534	–
904		me_i = atid_row(i)
905		me_j = atid_col(j)
537	–
906		#else
539	–
907		me_i = atid(i)
908		me_j = atid(j)
542	–
909		#endif
910	<
911	<	if (FF_uses_LJ .and. SimUsesLJ()) then
912	<	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
913	<	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
914	<
915	<	if ( is_LJ_i .and. is_LJ_j ) &
916	<	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
910	>
911	>	if (FF_uses_LJ .and. SIM_uses_LJ) then
912	>
913	>	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
914	>	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
915	>	endif
916	>
917		endif
918	<
919	<	if (FF_uses_dipoles .and. SimUsesDipoles()) then
554	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
555	<	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
918	>
919	>	if (FF_uses_charges .and. SIM_uses_charges) then
920
921	<	if ( is_DP_i .and. is_DP_j ) then
922	<
921	>	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
922	>	call do_charge_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
923	>	endif
924	>
925	>	endif
926	>
927	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
928	>
929	>	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
930		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
931		do_pot, do_stress)
932	<	if (FF_uses_RF .and. SimUsesRF()) then
932	>	if (FF_uses_RF .and. SIM_uses_RF) then
933		call accumulate_rf(i, j, r, u_l)
934		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
935	<	endif
565	<
935	>	endif
936		endif
937	+
938		endif
939
940	<	if (FF_uses_Sticky .and. SimUsesSticky()) then
940	>	if (FF_uses_Sticky .and. SIM_uses_sticky) then
941
942	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
572	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
573	<
574	<	if ( is_Sticky_i .and. is_Sticky_j ) then
942	>	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
943		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
944		do_pot, do_stress)
945		endif
946	+
947		endif
948
949
950	<	if (FF_uses_GB .and. SimUsesGB()) then
582	<
583	<	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
584	<	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
950	>	if (FF_uses_GB .and. SIM_uses_GB) then
951
952	<	if ( is_GB_i .and. is_GB_j ) then
952	>	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
953		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
954		do_pot, do_stress)
955		endif
956	+
957		endif
958	<
958	>
959	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
960	>
961	>	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
962	>	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
963	>	endif
964	>
965	>	endif
966	>
967		end subroutine do_pair
968
969
970	+
971	+	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
972	+	real( kind = dp ) :: pot
973	+	real( kind = dp ), dimension(3,nLocal) :: u_l
974	+	real (kind=dp), dimension(9,nLocal) :: A
975	+	real (kind=dp), dimension(3,nLocal) :: f
976	+	real (kind=dp), dimension(3,nLocal) :: t
977	+
978	+	logical, intent(inout) :: do_pot, do_stress
979	+	integer, intent(in) :: i, j
980	+	real ( kind = dp ), intent(inout) :: rijsq
981	+	real ( kind = dp ) :: r
982	+	real ( kind = dp ), intent(inout) :: d(3)
983	+
984	+	logical :: is_EAM_i, is_EAM_j
985	+
986	+	integer :: me_i, me_j
987	+
988	+	r = sqrt(rijsq)
989	+
990	+
991	+	#ifdef IS_MPI
992	+	if (tagRow(i) .eq. tagColumn(j)) then
993	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
994	+	endif
995	+
996	+	me_i = atid_row(i)
997	+	me_j = atid_col(j)
998	+
999	+	#else
1000	+
1001	+	me_i = atid(i)
1002	+	me_j = atid(j)
1003	+
1004	+	#endif
1005	+
1006	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
1007	+
1008	+	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1009	+	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1010	+
1011	+	endif
1012	+
1013	+	end subroutine do_prepair
1014	+
1015	+
1016	+
1017	+
1018	+	subroutine do_preforce(nlocal,pot)
1019	+	integer :: nlocal
1020	+	real( kind = dp ) :: pot
1021	+
1022	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
1023	+	call calc_EAM_preforce_Frho(nlocal,pot)
1024	+	endif
1025	+
1026	+
1027	+	end subroutine do_preforce
1028	+
1029	+
1030		subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1031
1032		real (kind = dp), dimension(3) :: q_i
1033		real (kind = dp), dimension(3) :: q_j
1034		real ( kind = dp ), intent(out) :: r_sq
1035	<	real( kind = dp ) :: d(3)
1036	<	real( kind = dp ) :: d_old(3)
1035	>	real( kind = dp ) :: d(3), scaled(3)
1036	>	integer i
1037	>
1038		d(1:3) = q_j(1:3) - q_i(1:3)
1039	<	d_old = d
1039	>
1040		! Wrap back into periodic box if necessary
1041	<	if ( SimUsesPBC() ) then
1041	>	if ( SIM_uses_PBC ) then
1042
1043	<	d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
1044	<	int(abs(d(1:3)/box(1:3)) + 0.5_dp)
1045	<
1046	<	endif
1047	<	r_sq = dot_product(d,d)
1048	<
613	<	end subroutine get_interatomic_vector
1043	>	if( .not.boxIsOrthorhombic ) then
1044	>	! calc the scaled coordinates.
1045	>
1046	>	scaled = matmul(HmatInv, d)
1047	>
1048	>	! wrap the scaled coordinates
1049
1050	<	subroutine check_initialization(error)
1051	<	integer, intent(out) :: error
617	<
618	<	error = 0
619	<	! Make sure we are properly initialized.
620	<	if (.not. do_forces_initialized) then
621	<	error = -1
622	<	return
623	<	endif
1050	>	scaled = scaled - anint(scaled)
1051	>
1052
1053	<	#ifdef IS_MPI
1054	<	if (.not. isMPISimSet()) then
1055	<	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
1056	<	error = -1
1057	<	return
1053	>	! calc the wrapped real coordinates from the wrapped scaled
1054	>	! coordinates
1055	>
1056	>	d = matmul(Hmat,scaled)
1057	>
1058	>	else
1059	>	! calc the scaled coordinates.
1060	>
1061	>	do i = 1, 3
1062	>	scaled(i) = d(i) * HmatInv(i,i)
1063	>
1064	>	! wrap the scaled coordinates
1065	>
1066	>	scaled(i) = scaled(i) - anint(scaled(i))
1067	>
1068	>	! calc the wrapped real coordinates from the wrapped scaled
1069	>	! coordinates
1070	>
1071	>	d(i) = scaled(i)*Hmat(i,i)
1072	>	enddo
1073	>	endif
1074	>
1075		endif
631	–	#endif
1076
1077	<	return
1078	<	end subroutine check_initialization
1079	<
1077	>	r_sq = dot_product(d,d)
1078	>
1079	>	end subroutine get_interatomic_vector
1080
1081		subroutine zero_work_arrays()
1082
#	Line 665 \| Line 1109 \| contains
1109
1110		#endif
1111
1112	+
1113	+	if (FF_uses_EAM .and. SIM_uses_EAM) then
1114	+	call clean_EAM()
1115	+	endif
1116	+
1117	+
1118	+
1119	+
1120	+
1121		rf = 0.0_dp
1122		tau_Temp = 0.0_dp
1123		virial_Temp = 0.0_dp
#	Line 777 \| Line 1230 \| end module do_Forces
1230		doesit = FF_uses_RF
1231		end function FF_RequiresPostpairCalc
1232
1233	+	#ifdef PROFILE
1234	+	function getforcetime() result(totalforcetime)
1235	+	real(kind=dp) :: totalforcetime
1236	+	totalforcetime = forcetime
1237	+	end function getforcetime
1238	+	#endif
1239	+
1240	+	!! This cleans componets of force arrays belonging only to fortran
1241	+
1242		end module do_Forces

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents): Revision 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC vs. Revision 1132 by tim, Sat Apr 24 04:31:36 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC vs.
Revision 1132 by tim, Sat Apr 24 04:31:36 2004 UTC