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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 441 by chuckv, Tue Apr 1 16:50:14 2003 UTC vs.
Revision 898 by chuckv, Mon Jan 5 22:49:14 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.7 2003-04-01 16:50:14 chuckv Exp $, $Date: 2003-04-01 16:50:14 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
7	>	!! @version $Id: do_Forces.F90,v 1.44 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
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 :: do_forces_initialized = .false., haveRlist = .false.
34	>	logical, save :: havePolicies = .false.
35		logical, save :: FF_uses_LJ
36		logical, save :: FF_uses_sticky
37		logical, save :: FF_uses_dipoles
#	Line 35 \| Line 39 \| module do_Forces
39		logical, save :: FF_uses_GB
40		logical, save :: FF_uses_EAM
41
42	+	real(kind=dp), save :: rlist, rlistsq
43	+
44		public :: init_FF
45		public :: do_force_loop
46	+	public :: setRlistDF
47
48	+	#ifdef PROFILE
49	+	public :: getforcetime
50	+	real, save :: forceTime = 0
51	+	real :: forceTimeInitial, forceTimeFinal
52	+	integer :: nLoops
53	+	#endif
54	+
55	+	logical, allocatable :: propertyMapI(:,:)
56	+	logical, allocatable :: propertyMapJ(:,:)
57	+
58		contains
59
60	+	subroutine setRlistDF( this_rlist )
61	+
62	+	real(kind=dp) :: this_rlist
63	+
64	+	rlist = this_rlist
65	+	rlistsq = rlist * rlist
66	+
67	+	haveRlist = .true.
68	+	if( havePolicies ) do_forces_initialized = .true.
69	+
70	+	end subroutine setRlistDF
71	+
72		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
73
74		integer, intent(in) :: LJMIXPOLICY
#	Line 87 \| Line 116 \| contains
116		!! check to make sure the FF_uses_RF setting makes sense
117
118		if (FF_uses_dipoles) then
90	–	rrf = getRrf()
91	–	rt = getRt()
92	–	call initialize_dipole(rrf, rt)
119		if (FF_uses_RF) then
120		dielect = getDielect()
121	<	call initialize_rf(rrf, rt, dielect)
121	>	call initialize_rf(dielect)
122		endif
123		else
124		if (FF_uses_RF) then
#	Line 100 \| Line 126 \| contains
126		thisStat = -1
127		return
128		endif
129	<	endif
129	>	endif
130
131		if (FF_uses_LJ) then
132
107	–	call getRcut(rcut)
108	–
133		select case (LJMIXPOLICY)
134		case (LB_MIXING_RULE)
135	<	call init_lj_FF(LB_MIXING_RULE, rcut, my_status)
135	>	call init_lj_FF(LB_MIXING_RULE, my_status)
136		case (EXPLICIT_MIXING_RULE)
137	<	call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
137	>	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
138		case default
139		write(default_error,*) 'unknown LJ Mixing Policy!'
140		thisStat = -1
#	Line 129 \| Line 153 \| contains
153		return
154		end if
155		endif
156	+
157	+
158	+	if (FF_uses_EAM) then
159	+	call init_EAM_FF(my_status)
160	+	if (my_status /= 0) then
161	+	write(,) "init_EAM_FF returned a bad status"
162	+	thisStat = -1
163	+	return
164	+	end if
165	+	endif
166	+
167	+
168
169		if (FF_uses_GB) then
170		call check_gb_pair_FF(my_status)
#	Line 140 \| Line 176 \| contains
176
177		if (FF_uses_GB .and. FF_uses_LJ) then
178		endif
179	+	if (.not. do_forces_initialized) then
180	+	!! Create neighbor lists
181	+	call expandNeighborList(nLocal, my_status)
182	+	if (my_Status /= 0) then
183	+	write(default_error,*) "SimSetup: ExpandNeighborList returned error."
184	+	thisStat = -1
185	+	return
186	+	endif
187	+	endif
188	+
189
190	+	havePolicies = .true.
191	+	if( haveRlist ) do_forces_initialized = .true.
192
145	–	do_forces_initialized = .true.
146	–
193		end subroutine init_FF
194
195
#	Line 152 \| Line 198 \| contains
198		subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
199		error)
200		!! Position array provided by C, dimensioned by getNlocal
201	<	real ( kind = dp ), dimension(3,getNlocal()) :: q
201	>	real ( kind = dp ), dimension(3,nLocal) :: q
202		!! Rotation Matrix for each long range particle in simulation.
203	<	real( kind = dp), dimension(9,getNlocal()) :: A
203	>	real( kind = dp), dimension(9,nLocal) :: A
204		!! Unit vectors for dipoles (lab frame)
205	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
205	>	real( kind = dp ), dimension(3,nLocal) :: u_l
206		!! Force array provided by C, dimensioned by getNlocal
207	<	real ( kind = dp ), dimension(3,getNlocal()) :: f
207	>	real ( kind = dp ), dimension(3,nLocal) :: f
208		!! Torsion array provided by C, dimensioned by getNlocal
209	<	real( kind = dp ), dimension(3,getNlocal()) :: t
209	>	real( kind = dp ), dimension(3,nLocal) :: t
210	>
211		!! Stress Tensor
212		real( kind = dp), dimension(9) :: tau
213		real ( kind = dp ) :: pot
#	Line 171 \| Line 218 \| contains
218		real( kind = DP ) :: pot_local
219		integer :: nrow
220		integer :: ncol
221	+	integer :: nprocs
222		#endif
175	–	integer :: nlocal
223		integer :: natoms
224		logical :: update_nlist
225		integer :: i, j, jbeg, jend, jnab
226		integer :: nlist
227	<	real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut
227	>	real( kind = DP ) :: rijsq
228		real(kind=dp),dimension(3) :: d
229		real(kind=dp) :: rfpot, mu_i, virial
230	<	integer :: me_i
230	>	integer :: me_i, me_j
231		logical :: is_dp_i
232		integer :: neighborListSize
233		integer :: listerror, error
234		integer :: localError
235	+	integer :: propPack_i, propPack_j
236
237	+	real(kind=dp) :: listSkin = 1.0
238	+
239		!! initialize local variables
240
241		#ifdef IS_MPI
242		pot_local = 0.0_dp
193	–	nlocal = getNlocal()
243		nrow = getNrow(plan_row)
244		ncol = getNcol(plan_col)
245		#else
197	–	nlocal = getNlocal()
246		natoms = nlocal
247		#endif
248	<
201	<	call getRcut(rcut,rc2=rcutsq)
202	<	call getRlist(rlist,rlistsq)
203	<
248	>
249		call check_initialization(localError)
250		if ( localError .ne. 0 ) then
251	+	call handleError("do_force_loop","Not Initialized")
252		error = -1
253		return
254		end if
#	Line 210 \| Line 256 \| contains
256
257		do_pot = do_pot_c
258		do_stress = do_stress_c
259	+
260	+
261	+	#ifdef IS_MPI
262	+	if (.not.allocated(propertyMapI)) then
263	+	allocate(propertyMapI(5,nrow))
264	+	endif
265	+
266	+	do i = 1, nrow
267	+	me_i = atid_row(i)
268	+	#else
269	+	if (.not.allocated(propertyMapI)) then
270	+	allocate(propertyMapI(5,nlocal))
271	+	endif
272	+
273	+	do i = 1, natoms
274	+	me_i = atid(i)
275	+	#endif
276	+
277	+	propertyMapI(1:5,i) = .false.
278	+
279	+	call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
280	+
281	+	! unpack the properties
282	+
283	+	if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
284	+	propertyMapI(1, i) = .true.
285	+	if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
286	+	propertyMapI(2, i) = .true.
287	+	if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
288	+	propertyMapI(3, i) = .true.
289	+	if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
290	+	propertyMapI(4, i) = .true.
291	+	if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
292	+	propertyMapI(5, i) = .true.
293	+
294	+	end do
295	+
296	+	#ifdef IS_MPI
297	+	if (.not.allocated(propertyMapJ)) then
298	+	allocate(propertyMapJ(5,ncol))
299	+	endif
300	+
301	+	do j = 1, ncol
302	+	me_j = atid_col(j)
303	+	#else
304	+	if (.not.allocated(propertyMapJ)) then
305	+	allocate(propertyMapJ(5,nlocal))
306	+	endif
307
308	+	do j = 1, natoms
309	+	me_j = atid(j)
310	+	#endif
311	+
312	+	propertyMapJ(1:5,j) = .false.
313	+
314	+	call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
315	+
316	+	! unpack the properties
317	+
318	+	if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
319	+	propertyMapJ(1, j) = .true.
320	+	if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
321	+	propertyMapJ(2, j) = .true.
322	+	if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
323	+	propertyMapJ(3, j) = .true.
324	+	if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
325	+	propertyMapJ(4, j) = .true.
326	+	if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
327	+	propertyMapJ(5, j) = .true.
328	+
329	+	end do
330	+
331		! Gather all information needed by all force loops:
332
333		#ifdef IS_MPI
#	Line 227 \| Line 344 \| contains
344		endif
345
346		#endif
347	<
347	>
348	>	!! Begin force loop timing:
349	>	#ifdef PROFILE
350	>	call cpu_time(forceTimeInitial)
351	>	nloops = nloops + 1
352	>	#endif
353	>
354		if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
355		!! See if we need to update neighbor lists
356	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
356	>	call checkNeighborList(nlocal, q, listSkin, update_nlist)
357		!! if_mpi_gather_stuff_for_prepair
358		!! do_prepair_loop_if_needed
359		!! if_mpi_scatter_stuff_from_prepair
360		!! if_mpi_gather_stuff_from_prepair_to_main_loop
361	<	else
362	<	!! See if we need to update neighbor lists
363	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
361	>
362	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
363	>	#ifdef IS_MPI
364	>
365	>	if (update_nlist) then
366	>
367	>	!! save current configuration, construct neighbor list,
368	>	!! and calculate forces
369	>	call saveNeighborList(nlocal, q)
370	>
371	>	neighborListSize = size(list)
372	>	nlist = 0
373	>
374	>	do i = 1, nrow
375	>	point(i) = nlist + 1
376	>
377	>	prepair_inner: do j = 1, ncol
378	>
379	>	if (skipThisPair(i,j)) cycle prepair_inner
380	>
381	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
382	>
383	>	if (rijsq < rlistsq) then
384	>
385	>	nlist = nlist + 1
386	>
387	>	if (nlist > neighborListSize) then
388	>	call expandNeighborList(nlocal, listerror)
389	>	if (listerror /= 0) then
390	>	error = -1
391	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
392	>	return
393	>	end if
394	>	neighborListSize = size(list)
395	>	endif
396	>
397	>	list(nlist) = j
398	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
399	>	endif
400	>	enddo prepair_inner
401	>	enddo
402	>
403	>	point(nrow + 1) = nlist + 1
404	>
405	>	else !! (of update_check)
406	>
407	>	! use the list to find the neighbors
408	>	do i = 1, nrow
409	>	JBEG = POINT(i)
410	>	JEND = POINT(i+1) - 1
411	>	! check thiat molecule i has neighbors
412	>	if (jbeg .le. jend) then
413	>
414	>	do jnab = jbeg, jend
415	>	j = list(jnab)
416	>
417	>	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
418	>	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
419	>	u_l, A, f, t, pot_local)
420	>
421	>	enddo
422	>	endif
423	>	enddo
424		endif
425
426	<	#ifdef IS_MPI
426	>	#else
427
428		if (update_nlist) then
429
430	+	! save current configuration, contruct neighbor list,
431	+	! and calculate forces
432	+	call saveNeighborList(natoms, q)
433	+
434	+	neighborListSize = size(list)
435	+
436	+	nlist = 0
437	+
438	+	do i = 1, natoms-1
439	+	point(i) = nlist + 1
440	+
441	+	prepair_inner: do j = i+1, natoms
442	+
443	+	if (skipThisPair(i,j)) cycle prepair_inner
444	+
445	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
446	+
447	+
448	+	if (rijsq < rlistsq) then
449	+
450	+
451	+	nlist = nlist + 1
452	+
453	+	if (nlist > neighborListSize) then
454	+	call expandNeighborList(natoms, listerror)
455	+	if (listerror /= 0) then
456	+	error = -1
457	+	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
458	+	return
459	+	end if
460	+	neighborListSize = size(list)
461	+	endif
462	+
463	+	list(nlist) = j
464	+
465	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
466	+	u_l, A, f, t, pot)
467	+
468	+	endif
469	+	enddo prepair_inner
470	+	enddo
471	+
472	+	point(natoms) = nlist + 1
473	+
474	+	else !! (update)
475	+
476	+	! use the list to find the neighbors
477	+	do i = 1, natoms-1
478	+	JBEG = POINT(i)
479	+	JEND = POINT(i+1) - 1
480	+	! check thiat molecule i has neighbors
481	+	if (jbeg .le. jend) then
482	+
483	+	do jnab = jbeg, jend
484	+	j = list(jnab)
485	+
486	+	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
487	+	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
488	+	u_l, A, f, t, pot)
489	+
490	+	enddo
491	+	endif
492	+	enddo
493	+	endif
494	+	#endif
495	+	!! Do rest of preforce calculations
496	+	!! do necessary preforce calculations
497	+	call do_preforce(nlocal,pot)
498	+	! we have already updated the neighbor list set it to false...
499	+	update_nlist = .false.
500	+	else
501	+	!! See if we need to update neighbor lists for non pre-pair
502	+	call checkNeighborList(nlocal, q, listSkin, update_nlist)
503	+	endif
504	+
505	+
506	+
507	+
508	+
509	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
510	+
511	+
512	+
513	+
514	+
515	+	#ifdef IS_MPI
516	+
517	+	if (update_nlist) then
518		!! save current configuration, construct neighbor list,
519		!! and calculate forces
520	<	call saveNeighborList(q)
520	>	call saveNeighborList(nlocal, q)
521
522		neighborListSize = size(list)
523		nlist = 0
524
525		do i = 1, nrow
526	+
527		point(i) = nlist + 1
528
529		inner: do j = 1, ncol
#	Line 260 \| Line 532 \| contains
532
533		call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
534
535	<	if (rijsq < rlistsq) then
535	>	if (rijsq < rlistsq) then
536
537		nlist = nlist + 1
538
#	Line 276 \| Line 548 \| contains
548
549		list(nlist) = j
550
551	<	if (rijsq < rcutsq) then
552	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
553	<	u_l, A, f, t, pot_local)
282	<	endif
551	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
552	>	u_l, A, f, t, pot_local)
553	>
554		endif
555		enddo inner
556		enddo
#	Line 310 \| Line 581 \| contains
581		#else
582
583		if (update_nlist) then
584	<
584	>
585		! save current configuration, contruct neighbor list,
586		! and calculate forces
587	<	call saveNeighborList(q)
587	>	call saveNeighborList(natoms, q)
588
589		neighborListSize = size(list)
590
#	Line 329 \| Line 600 \| contains
600		call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
601
602
603	<	if (rijsq < rlistsq) then
603	>	if (rijsq < rlistsq) then
604
605		nlist = nlist + 1
606
#	Line 345 \| Line 616 \| contains
616
617		list(nlist) = j
618
619	<	if (rijsq < rcutsq) then
349	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
619	>	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
620		u_l, A, f, t, pot)
621	<	endif
621	>
622		endif
623		enddo inner
624		enddo
#	Line 379 \| Line 649 \| contains
649		#endif
650
651		! phew, done with main loop.
652	<
652	>
653	>	!! Do timing
654	>	#ifdef PROFILE
655	>	call cpu_time(forceTimeFinal)
656	>	forceTime = forceTime + forceTimeFinal - forceTimeInitial
657	>	#endif
658	>
659	>
660		#ifdef IS_MPI
661		!!distribute forces
662
#	Line 441 \| Line 718 \| contains
718		end do
719		#endif
720
721	<	do i = 1, getNlocal()
721	>	do i = 1, nLocal
722
723		rfpot = 0.0_DP
724		#ifdef IS_MPI
#	Line 479 \| Line 756 \| contains
756		endif
757
758		if (do_stress) then
759	<	call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
759	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
760		mpi_comm_world,mpi_err)
761	<	call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
761	>	call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
762		mpi_comm_world,mpi_err)
763		endif
764
#	Line 491 \| Line 768 \| contains
768		tau = tau_Temp
769		virial = virial_Temp
770		endif
771	<
771	>
772		#endif
773
774	+
775	+
776		end subroutine do_force_loop
777
778		subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
779
780		real( kind = dp ) :: pot
781	<	real( kind = dp ), dimension(:,:) :: u_l
782	<	real (kind=dp), dimension(:,:) :: A
783	<	real (kind=dp), dimension(:,:) :: f
784	<	real (kind=dp), dimension(:,:) :: t
781	>	real( kind = dp ), dimension(3,nLocal) :: u_l
782	>	real (kind=dp), dimension(9,nLocal) :: A
783	>	real (kind=dp), dimension(3,nLocal) :: f
784	>	real (kind=dp), dimension(3,nLocal) :: t
785
786		logical, intent(inout) :: do_pot, do_stress
787		integer, intent(in) :: i, j
#	Line 512 \| Line 791 \| contains
791		logical :: is_LJ_i, is_LJ_j
792		logical :: is_DP_i, is_DP_j
793		logical :: is_GB_i, is_GB_j
794	+	logical :: is_EAM_i,is_EAM_j
795		logical :: is_Sticky_i, is_Sticky_j
796		integer :: me_i, me_j
797	<
797	>	integer :: propPack_i
798	>	integer :: propPack_j
799		r = sqrt(rijsq)
800
801		#ifdef IS_MPI
802	+	if (tagRow(i) .eq. tagColumn(j)) then
803	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
804	+	endif
805
806		me_i = atid_row(i)
807		me_j = atid_col(j)
#	Line 528 \| Line 812 \| contains
812		me_j = atid(j)
813
814		#endif
815	<
815	>
816		if (FF_uses_LJ .and. SimUsesLJ()) then
533	–	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
534	–	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
817
818	<	if ( is_LJ_i .and. is_LJ_j ) &
818	>	if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
819		call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
820	+
821		endif
822
823		if (FF_uses_dipoles .and. SimUsesDipoles()) then
541	–	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
542	–	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
824
825	<	if ( is_DP_i .and. is_DP_j ) then
545	<
825	>	if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
826		call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
827		do_pot, do_stress)
828		if (FF_uses_RF .and. SimUsesRF()) then
#	Line 555 \| Line 835 \| contains
835
836		if (FF_uses_Sticky .and. SimUsesSticky()) then
837
838	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
559	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
560	<
561	<	if ( is_Sticky_i .and. is_Sticky_j ) then
838	>	if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
839		call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
840		do_pot, do_stress)
841		endif
#	Line 566 \| Line 843 \| contains
843
844
845		if (FF_uses_GB .and. SimUsesGB()) then
569	–
570	–	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
571	–	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
846
847	<	if ( is_GB_i .and. is_GB_j ) then
847	>	if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
848		call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
849		do_pot, do_stress)
850		endif
851	+
852		endif
853
854	+
855	+
856	+	if (FF_uses_EAM .and. SimUsesEAM()) then
857	+
858	+	if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
859	+	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
860	+	endif
861	+
862	+	endif
863	+
864		end subroutine do_pair
865
866
867	+
868	+	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
869	+	real( kind = dp ) :: pot
870	+	real( kind = dp ), dimension(3,nLocal) :: u_l
871	+	real (kind=dp), dimension(9,nLocal) :: A
872	+	real (kind=dp), dimension(3,nLocal) :: f
873	+	real (kind=dp), dimension(3,nLocal) :: t
874	+
875	+	logical, intent(inout) :: do_pot, do_stress
876	+	integer, intent(in) :: i, j
877	+	real ( kind = dp ), intent(inout) :: rijsq
878	+	real ( kind = dp ) :: r
879	+	real ( kind = dp ), intent(inout) :: d(3)
880	+
881	+	logical :: is_EAM_i, is_EAM_j
882	+
883	+	integer :: me_i, me_j
884	+
885	+	r = sqrt(rijsq)
886	+
887	+
888	+	#ifdef IS_MPI
889	+	if (tagRow(i) .eq. tagColumn(j)) then
890	+	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
891	+	endif
892	+
893	+	me_i = atid_row(i)
894	+	me_j = atid_col(j)
895	+
896	+	#else
897	+
898	+	me_i = atid(i)
899	+	me_j = atid(j)
900	+
901	+	#endif
902	+
903	+	if (FF_uses_EAM .and. SimUsesEAM()) then
904	+	call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
905	+	call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
906	+
907	+	if ( is_EAM_i .and. is_EAM_j ) &
908	+	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
909	+	endif
910	+
911	+	end subroutine do_prepair
912	+
913	+
914	+
915	+
916	+	subroutine do_preforce(nlocal,pot)
917	+	integer :: nlocal
918	+	real( kind = dp ) :: pot
919	+
920	+	if (FF_uses_EAM .and. SimUsesEAM()) then
921	+	call calc_EAM_preforce_Frho(nlocal,pot)
922	+	endif
923	+
924	+
925	+	end subroutine do_preforce
926	+
927	+
928		subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
929
930		real (kind = dp), dimension(3) :: q_i
931		real (kind = dp), dimension(3) :: q_j
932		real ( kind = dp ), intent(out) :: r_sq
933	<	real( kind = dp ) :: d(3)
934	<	real( kind = dp ) :: d_old(3)
935	<	d(1:3) = q_i(1:3) - q_j(1:3)
936	<	d_old = d
933	>	real( kind = dp ) :: d(3), scaled(3)
934	>	integer i
935	>
936	>	d(1:3) = q_j(1:3) - q_i(1:3)
937	>
938		! Wrap back into periodic box if necessary
939		if ( SimUsesPBC() ) then
940
941	<	d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
942	<	int(abs(d(1:3)/box(1:3)) + 0.5_dp)
941	>	if( .not.boxIsOrthorhombic ) then
942	>	! calc the scaled coordinates.
943	>
944	>	scaled = matmul(HmatInv, d)
945	>
946	>	! wrap the scaled coordinates
947	>
948	>	scaled = scaled - anint(scaled)
949	>
950	>
951	>	! calc the wrapped real coordinates from the wrapped scaled
952	>	! coordinates
953	>
954	>	d = matmul(Hmat,scaled)
955	>
956	>	else
957	>	! calc the scaled coordinates.
958	>
959	>	do i = 1, 3
960	>	scaled(i) = d(i) * HmatInv(i,i)
961	>
962	>	! wrap the scaled coordinates
963	>
964	>	scaled(i) = scaled(i) - anint(scaled(i))
965	>
966	>	! calc the wrapped real coordinates from the wrapped scaled
967	>	! coordinates
968	>
969	>	d(i) = scaled(i)*Hmat(i,i)
970	>	enddo
971	>	endif
972
973		endif
974	+
975		r_sq = dot_product(d,d)
976	<
976	>
977		end subroutine get_interatomic_vector
978	<
978	>
979		subroutine check_initialization(error)
980		integer, intent(out) :: error
981
982		error = 0
983		! Make sure we are properly initialized.
984		if (.not. do_forces_initialized) then
985	+	write(,) "Forces not initialized"
986		error = -1
987		return
988		endif
#	Line 652 \| Line 1030 \| contains
1030
1031		#endif
1032
1033	+
1034	+	if (FF_uses_EAM .and. SimUsesEAM()) then
1035	+	call clean_EAM()
1036	+	endif
1037	+
1038	+
1039	+
1040	+
1041	+
1042		rf = 0.0_dp
1043		tau_Temp = 0.0_dp
1044		virial_Temp = 0.0_dp
658	–
1045		end subroutine zero_work_arrays
1046
1047		function skipThisPair(atom1, atom2) result(skip_it)
#	Line 765 \| Line 1151 \| end module do_Forces
1151		doesit = FF_uses_RF
1152		end function FF_RequiresPostpairCalc
1153
1154	+	#ifdef PROFILE
1155	+	function getforcetime() result(totalforcetime)
1156	+	real(kind=dp) :: totalforcetime
1157	+	totalforcetime = forcetime
1158	+	end function getforcetime
1159	+	#endif
1160	+
1161	+	!! This cleans componets of force arrays belonging only to fortran
1162	+
1163		end module do_Forces

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents): Revision 441 by chuckv, Tue Apr 1 16:50:14 2003 UTC vs. Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 441 by chuckv, Tue Apr 1 16:50:14 2003 UTC vs.
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC