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 1113 by tim, Thu Apr 15 16:18:26 2004 UTC vs.
Revision 1150 by gezelter, Fri May 7 21:35:05 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.49 2004-04-15 16:18:26 tim Exp $, $Date: 2004-04-15 16:18:26 $, $Name: not supported by cvs2svn $, $Revision: 1.49 $
7	>	!! @version $Id: do_Forces.F90,v 1.54 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
8
9		module do_Forces
10		use force_globals
11		use simulation
12		use definitions
13		use atype_module
14	+	use switcheroo
15		use neighborLists
16		use lj
17		use sticky_pair
#	Line 30 | Line 31 | module do_Forces
31
32		#define __FORTRAN90
33		#include "fForceField.h"
34	+	#include "fSwitchingFunction.h"
35
36		logical, save :: haveRlist = .false.
37		logical, save :: haveNeighborList = .false.
#	Line 55 | Line 57 | module do_Forces
57		logical, save :: SIM_requires_prepair_calc
58		logical, save :: SIM_uses_directional_atoms
59		logical, save :: SIM_uses_PBC
60	+	logical, save :: SIM_uses_molecular_cutoffs
61
62		real(kind=dp), save :: rlist, rlistsq
63
64		public :: init_FF
65		public :: do_force_loop
66		public :: setRlistDF
64	–
67
68		#ifdef PROFILE
69		public :: getforcetime
#	Line 160 | Line 162 | contains
162		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
163		SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
164		SIM_uses_PBC = SimUsesPBC()
165	+	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
166
167		haveSIMvariables = .true.
168
#	Line 277 | Line 280 | contains
280
281		!! Assume sanity (for the sake of argument)
282		haveSaneForceField = .true.
283	+	!!
284	+	if (FF_uses_charges) then
285	+	dielect = getDielect()
286	+	call initialize_charge(dielect)
287	+	endif
288	+
289
290		!! check to make sure the FF_uses_RF setting makes sense
291
#	Line 356 | Line 365 | contains
365		endif
366		haveNeighborList = .true.
367		endif
368	+
369
370	+
371		end subroutine init_FF
372
373
374		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
375		!------------------------------------------------------------->
376	<	subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
377	<	error)
376	>	subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
377	>	do_pot_c, do_stress_c, error)
378		!! Position array provided by C, dimensioned by getNlocal
379		real ( kind = dp ), dimension(3,nLocal) :: q
380	+	!! molecular center-of-mass position array
381	+	real ( kind = dp ), dimension(3,nGroup) :: q_group
382		!! Rotation Matrix for each long range particle in simulation.
383		real( kind = dp), dimension(9,nLocal) :: A
384		!! Unit vectors for dipoles (lab frame)
#	Line 381 | Line 394 | contains
394		logical ( kind = 2) :: do_pot_c, do_stress_c
395		logical :: do_pot
396		logical :: do_stress
397	+	logical :: in_switching_region
398		#ifdef IS_MPI
399		real( kind = DP ) :: pot_local
400		integer :: nrow
401		integer :: ncol
402		integer :: nprocs
403	+	integer :: nrow_group
404	+	integer :: ncol_group
405		#endif
406		integer :: natoms
407		logical :: update_nlist
408		integer :: i, j, jbeg, jend, jnab
409	+	integer :: ia, jb, atom1, atom2
410		integer :: nlist
411	<	real( kind = DP ) ::  rijsq
412	<	real(kind=dp),dimension(3) :: d
411	>	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vab
412	>	real( kind = DP ) :: sw, dswdr, swderiv, mf
413	>	real(kind=dp),dimension(3) :: d_atm, d_grp
414		real(kind=dp) :: rfpot, mu_i, virial
415		integer :: me_i, me_j
416		logical :: is_dp_i
#	Line 402 | Line 420 | contains
420		integer :: propPack_i, propPack_j
421
422		real(kind=dp) :: listSkin = 1.0
423	<
423	>
424		!! initialize local variables
425	<
425	>
426		#ifdef IS_MPI
427		pot_local = 0.0_dp
428		nrow   = getNrow(plan_row)
429		ncol   = getNcol(plan_col)
430	+	nrow_group   = getNrowGroup(plan_row)
431	+	ncol_group   = getNcolGroup(plan_col)
432		#else
433		natoms = nlocal
434		#endif
435	<
435	>
436		call doReadyCheck(localError)
437		if ( localError .ne. 0 ) then
438		call handleError("do_force_loop", "Not Initialized")
#	Line 420 | Line 440 | contains
440		return
441		end if
442		call zero_work_arrays()
443	<
443	>
444		do_pot = do_pot_c
445		do_stress = do_stress_c
446	<
446	>
447		! Gather all information needed by all force loops:
448
449		#ifdef IS_MPI
450	+
451	+	call gather(q, q_Row, plan_row3d)
452	+	call gather(q, q_Col, plan_col3d)
453
454	<	call gather(q,q_Row,plan_row3d)
455	<	call gather(q,q_Col,plan_col3d)
454	>	call gather(q_group, q_group_Row, plan_row_Group_3d)
455	>	call gather(q_group, q_group_Col, plan_col_Group_3d)
456
457		if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
458		call gather(u_l,u_l_Row,plan_row3d)
#	Line 440 | Line 463 | contains
463		endif
464
465		#endif
466	<
467	<	!! Begin force loop timing:
466	>
467	>	!! Begin force loop timing:
468		#ifdef PROFILE
469		call cpu_time(forceTimeInitial)
470		nloops = nloops + 1
471		#endif
472	<
472	>
473		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
474		!! See if we need to update neighbor lists
475	<	call checkNeighborList(nlocal, q, listSkin, update_nlist)
475	>
476	>	call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
477	>
478		!! if_mpi_gather_stuff_for_prepair
479		!! do_prepair_loop_if_needed
480		!! if_mpi_scatter_stuff_from_prepair
481		!! if_mpi_gather_stuff_from_prepair_to_main_loop
482	<
483	<	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
482	>
483	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
484		#ifdef IS_MPI
460	–
461	–	if (update_nlist) then
485
486	<	!! save current configuration, construct neighbor list,
464	<	!! and calculate forces
465	<	call saveNeighborList(nlocal, q)
466	<
467	<	neighborListSize = size(list)
468	<	nlist = 0
469	<
470	<	do i = 1, nrow
471	<	point(i) = nlist + 1
486	>	if (update_nlist) then
487
488	<	prepair_inner: do j = 1, ncol
489	<
490	<	if (skipThisPair(i,j)) cycle prepair_inner
491	<
492	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
493	<
494	<	if (rijsq < rlistsq) then
495	<
496	<	nlist = nlist + 1
488	>	!! save current configuration, construct neighbor list,
489	>	!! and calculate forces
490	>
491	>	call saveNeighborList(nGroup, q_group)
492	>
493	>	neighborListSize = size(list)
494	>	nlist = 0
495	>
496	>	do i = 1, nrow_group
497	>	point(i) = nlist + 1
498	>
499	>	do j = 1, ncol_group
500
501	<	if (nlist > neighborListSize) then
502	<	call expandNeighborList(nlocal, listerror)
485	<	if (listerror /= 0) then
486	<	error = -1
487	<	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
488	<	return
489	<	end if
490	<	neighborListSize = size(list)
491	<	endif
501	>	call get_interatomic_vector(q_group_Row(:,i), &
502	>	q_group_Col(:,j), d_grp, rgrpsq)
503
504	<	list(nlist) = j
505	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)
506	<	endif
507	<	enddo prepair_inner
508	<	enddo
504	>	if (rgrpsq < rlistsq) then
505	>	nlist = nlist + 1
506	>
507	>	if (nlist > neighborListSize) then
508	>	call expandNeighborList(nGroup, listerror)
509	>	if (listerror /= 0) then
510	>	error = -1
511	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
512	>	return
513	>	end if
514	>	neighborListSize = size(list)
515	>	endif
516	>
517	>	list(nlist) = j
518	>
519	>	do ia = groupStart(i), groupStart(i+1)-1
520	>	atom1 = groupList(ia)
521
522	<	point(nrow + 1) = nlist + 1
523	<
524	<	else  !! (of update_check)
522	>	prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
523	>	atom2 = groupList(jb)
524	>
525	>	if (skipThisPair(atom1, atom2)) cycle prepair_inner
526	>
527	>	call get_interatomic_vector(q_Row(:,atom1), &
528	>	q_Col(:,atom2), d_atm, ratmsq)
529	>
530	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
531	>	rgrpsq, d_grp, do_pot, do_stress, &
532	>	u_l, A, f, t, pot_local)
533
534	<	! use the list to find the neighbors
535	<	do i = 1, nrow
536	<	JBEG = POINT(i)
506	<	JEND = POINT(i+1) - 1
507	<	! check thiat molecule i has neighbors
508	<	if (jbeg .le. jend) then
509	<
510	<	do jnab = jbeg, jend
511	<	j = list(jnab)
512	<
513	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
514	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
515	<	u_l, A, f, t, pot_local)
516	<
534	>	enddo prepair_inner
535	>	enddo
536	>	end if
537		enddo
538	<	endif
539	<	enddo
540	<	endif
541	<
542	<	#else
543	<
544	<	if (update_nlist) then
538	>	enddo
539	>	point(nrow_group + 1) = nlist + 1
540	>
541	>	else  !! (of update_check)
542	>
543	>	! use the list to find the neighbors
544	>	do i = 1, nrow_group
545	>	JBEG = POINT(i)
546	>	JEND = POINT(i+1) - 1
547	>	! check that group i has neighbors
548	>	if (jbeg .le. jend) then
549	>
550	>	do jnab = jbeg, jend
551	>	j = list(jnab)
552	>
553	>	do ia = groupStart(i), groupStart(i+1)-1
554	>	atom1 = groupList(ia)
555	>
556	>	do jb = groupStart(j), groupStart(j+1)-1
557	>	atom2 = groupList(jb)
558	>
559	>	call get_interatomic_vector(q_Row(:,atom1), &
560	>	q_Col(:,atom2), d_atm, ratmsq)
561	>
562	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
563	>	rgrpsq, d_grp, do_pot, do_stress, &
564	>	u_l, A, f, t, pot_local)
565	>
566	>	enddo
567	>	enddo
568	>	enddo
569	>	endif
570	>	enddo
571	>	endif
572
573	<	! save current configuration, contruct neighbor list,
527	<	! and calculate forces
528	<	call saveNeighborList(natoms, q)
529	<
530	<	neighborListSize = size(list)
531	<
532	<	nlist = 0
573	>	#else
574
575	<	do i = 1, natoms-1
535	<	point(i) = nlist + 1
575	>	if (update_nlist) then
576
577	<	prepair_inner: do j = i+1, natoms
577	>	!! save current configuration, construct neighbor list,
578	>	!! and calculate forces
579	>
580	>	call saveNeighborList(nGroup, q_group)
581	>
582	>	neighborListSize = size(list)
583	>	nlist = 0
584	>
585	>	do i = 1, nGroup-1
586	>	point(i) = nlist + 1
587
588	<	if (skipThisPair(i,j))  cycle prepair_inner
540	<
541	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
542	<
543	<
544	<	if (rijsq < rlistsq) then
545	<
546	<
547	<	nlist = nlist + 1
548	<
549	<	if (nlist > neighborListSize) then
550	<	call expandNeighborList(natoms, listerror)
551	<	if (listerror /= 0) then
552	<	error = -1
553	<	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
554	<	return
555	<	end if
556	<	neighborListSize = size(list)
557	<	endif
588	>	do j = i+1, nGroup
589
590	<	list(nlist) = j
590	>	call get_interatomic_vector(q_group(:,i), &
591	>	q_group(:,j), d_grp, rgrpsq)
592
593	<	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
594	<	u_l, A, f, t, pot)
593	>	if (rgrpsq < rlistsq) then
594	>	nlist = nlist + 1
595	>
596	>	if (nlist > neighborListSize) then
597	>	call expandNeighborList(nGroup, listerror)
598	>	if (listerror /= 0) then
599	>	error = -1
600	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
601	>	return
602	>	end if
603	>	neighborListSize = size(list)
604	>	endif
605	>
606	>	list(nlist) = j
607	>
608	>	do ia = groupStart(i), groupStart(i+1)-1
609	>	atom1 = groupList(ia)
610	>
611	>	prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
612	>	atom2 = groupList(jb)
613	>
614	>	if (skipThisPair(atom1, atom2)) cycle prepair_inner
615	>
616	>	call get_interatomic_vector(q(:,atom1), &
617	>	q(:,atom2), d_atm, ratmsq)
618	>
619	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
620	>	rgrpsq, d_grp, do_pot, do_stress, &
621	>	u_l, A, f, t, pot)
622	>
623	>	enddo prepair_inner
624	>	enddo
625	>	end if
626	>	enddo
627	>	enddo
628	>	point(nGroup) = nlist + 1
629	>
630	>	else  !! (of update_check)
631	>
632	>	! use the list to find the neighbors
633	>	do i = 1, nGroup-1
634	>	JBEG = POINT(i)
635	>	JEND = POINT(i+1) - 1
636	>	! check that group i has neighbors
637	>	if (jbeg .le. jend) then
638
639	+	do jnab = jbeg, jend
640	+	j = list(jnab)
641	+
642	+	do ia = groupStart(i), groupStart(i+1)-1
643	+	atom1 = groupList(ia)
644	+
645	+	do jb = groupStart(j), groupStart(j+1)-1
646	+	atom2 = groupList(jb)
647	+
648	+	call get_interatomic_vector(q(:,atom1), &
649	+	q(:,atom2), d_atm, ratmsq)
650	+
651	+	call do_prepair(atom1, atom2, ratmsq, d_atm, &
652	+	rgrpsq, d_grp, do_pot, do_stress, &
653	+	u_l, A, f, t, pot)
654	+
655	+	enddo
656	+	enddo
657	+	enddo
658		endif
659	<	enddo prepair_inner
660	<	enddo
567	<
568	<	point(natoms) = nlist + 1
659	>	enddo
660	>	endif
661
570	–	else !! (update)
571	–
572	–	! use the list to find the neighbors
573	–	do i = 1, natoms-1
574	–	JBEG = POINT(i)
575	–	JEND = POINT(i+1) - 1
576	–	! check thiat molecule i has neighbors
577	–	if (jbeg .le. jend) then
578	–
579	–	do jnab = jbeg, jend
580	–	j = list(jnab)
581	–
582	–	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
583	–	call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
584	–	u_l, A, f, t, pot)
585	–
586	–	enddo
587	–	endif
588	–	enddo
589	–	endif
662		#endif
663	<	!! Do rest of preforce calculations
664	<	!! do necessary preforce calculations
665	<	call do_preforce(nlocal,pot)
666	<	! we have already updated the neighbor list set it to false...
667	<	update_nlist = .false.
663	>
664	>	!! Do rest of preforce calculations
665	>	!! do necessary preforce calculations
666	>	call do_preforce(nlocal,pot)
667	>	! we have already updated the neighbor list set it to false...
668	>	update_nlist = .false.
669		else
670		!! See if we need to update neighbor lists for non pre-pair
671	<	call checkNeighborList(nlocal, q, listSkin, update_nlist)
671	>	call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
672		endif
673	<
674	<
675	<
603	<
604	<
605	<	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
606	<
607	<
608	<
609	<
610	<
673	>
674	>	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
675	>
676		#ifdef IS_MPI
677
678		if (update_nlist) then
679	+
680		!! save current configuration, construct neighbor list,
681		!! and calculate forces
616	–	call saveNeighborList(nlocal, q)
682
683	+	call saveNeighborList(nGroup, q_group)
684	+
685		neighborListSize = size(list)
686		nlist = 0
687
688	<	do i = 1, nrow
622	<
688	>	do i = 1, nrow_group
689		point(i) = nlist + 1
690
691	<	inner: do j = 1, ncol
691	>	do j = 1, ncol_group
692
693	<	if (skipThisPair(i,j)) cycle inner
693	>	call get_interatomic_vector(q_group_Row(:,i), &
694	>	q_group_Col(:,j), d_grp, rgrpsq)
695
696	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
630	<
631	<	if (rijsq < rlistsq) then
632	<
696	>	if (rgrpsq < rlistsq) then
697		nlist = nlist + 1
698
699		if (nlist > neighborListSize) then
700	<	call expandNeighborList(nlocal, listerror)
700	>	call expandNeighborList(nGroup, listerror)
701		if (listerror /= 0) then
702		error = -1
703		write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
#	Line 643 | Line 707 | contains
707		endif
708
709		list(nlist) = j
646	–
647	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
648	–	u_l, A, f, t, pot_local)
710
711	<	endif
712	<	enddo inner
713	<	enddo
711	>	vab = 0.0d0
712	>	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
713	>	in_switching_region)
714	>
715	>	do ia = groupStart(i), groupStart(i+1)-1
716	>	atom1 = groupList(ia)
717	>
718	>	inner: do jb = groupStart(j), groupStart(j+1)-1
719	>	atom2 = groupList(jb)
720	>
721	>	if (skipThisPair(atom1, atom2)) cycle inner
722	>
723	>	call get_interatomic_vector(q_Row(:,atom1), &
724	>	q_Col(:,atom2), d_atm, ratmsq)
725	>
726	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
727	>	do_pot, do_stress, &
728	>	u_l, A, f, t, pot_local, vpair)
729	>
730	>	vab = vab + vpair
731	>	enddo inner
732	>	enddo
733
734	<	point(nrow + 1) = nlist + 1
734	>	if (in_switching_region) then
735	>	swderiv = vab*dswdr/rgrp
736	>
737	>	do ia=groupStart(i), groupStart(i+1)-1
738	>	atom1=groupList(ia)
739	>	mf = mfact(atom1)
740	>	f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
741	>	f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
742	>	f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
743	>	enddo
744	>
745	>	do jb=groupStart(j), groupStart(j+1)-1
746	>	atom2=groupList(jb)
747	>	mf = mfact(atom2)
748	>	f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
749	>	f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
750	>	f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
751	>	enddo
752	>	endif
753	>
754	>	end if
755	>	enddo
756	>	enddo
757	>	point(nrow_group + 1) = nlist + 1
758
759		else  !! (of update_check)
760	<
760	>
761		! use the list to find the neighbors
762	<	do i = 1, nrow
762	>	do i = 1, nrow_group
763		JBEG = POINT(i)
764		JEND = POINT(i+1) - 1
765	<	! check thiat molecule i has neighbors
765	>	! check that group i has neighbors
766		if (jbeg .le. jend) then
767
768		do jnab = jbeg, jend
769		j = list(jnab)
770	+
771	+	call get_interatomic_vector(q_group_Row(:,i), &
772	+	q_group_Col(:,j), d_grp, rgrpsq)
773	+
774	+	vab = 0.0d0
775	+	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
776	+	in_switching_region)
777	+
778	+	do ia = groupStart(i), groupStart(i+1)-1
779	+	atom1 = groupList(ia)
780	+
781	+	do jb = groupStart(j), groupStart(j+1)-1
782	+	atom2 = groupList(jb)
783	+
784	+	call get_interatomic_vector(q_Row(:,atom1), &
785	+	q_Col(:,atom2), d_atm, ratmsq)
786	+
787	+	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
788	+	do_pot, do_stress, &
789	+	u_l, A, f, t, pot_local, vpair)
790	+
791	+	vab = vab + vpair
792	+
793	+	enddo
794	+	enddo
795	+
796	+	if (in_switching_region) then
797	+	swderiv = vab*dswdr/rgrp
798	+
799	+	do ia=groupStart(i), groupStart(i+1)-1
800	+	atom1=groupList(ia)
801	+	mf = mfact(atom1)
802	+	f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
803	+	f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
804	+	f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
805	+	enddo
806	+
807	+	do jb=groupStart(j), groupStart(j+1)-1
808	+	atom2=groupList(jb)
809	+	mf = mfact(atom2)
810	+	f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
811	+	f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
812	+	f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
813	+	enddo
814	+	endif
815
668	–	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
669	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
670	–	u_l, A, f, t, pot_local)
671	–
816		enddo
817		endif
818		enddo
819		endif
820
821		#else
678	–
679	–	if (update_nlist) then
822
823	<	! save current configuration, contruct neighbor list,
682	<	! and calculate forces
683	<	call saveNeighborList(natoms, q)
823	>	if (update_nlist) then
824
825	+	!! save current configuration, construct neighbor list,
826	+	!! and calculate forces
827	+
828	+	call saveNeighborList(nGroup, q_group)
829	+
830		neighborListSize = size(list)
831	<
687	<	nlist = 0
831	>	nlist = 0
832
833	<	do i = 1, natoms-1
833	>	do i = 1, nGroup-1
834		point(i) = nlist + 1
835
836	<	inner: do j = i+1, natoms
836	>	do j = i+1, nGroup
837
838	<	if (skipThisPair(i,j))  cycle inner
839	<
840	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
841	<
698	<
699	<	if (rijsq < rlistsq) then
700	<
838	>	call get_interatomic_vector(q_group(:,i), &
839	>	q_group(:,j), d_grp, rgrpsq)
840	>
841	>	if (rgrpsq < rlistsq) then
842		nlist = nlist + 1
843	<
843	>
844		if (nlist > neighborListSize) then
845	<	call expandNeighborList(natoms, listerror)
845	>	call expandNeighborList(nGroup, listerror)
846		if (listerror /= 0) then
847		error = -1
848		write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
#	Line 711 | Line 852 | contains
852		endif
853
854		list(nlist) = j
855	+
856	+	vab = 0.0d0
857	+	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
858	+	in_switching_region)
859
860	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
861	<	u_l, A, f, t, pot)
862	<
863	<	endif
864	<	enddo inner
860	>	do ia = groupStart(i), groupStart(i+1)-1
861	>	atom1 = groupList(ia)
862	>
863	>	inner: do jb = groupStart(j), groupStart(j+1)-1
864	>	atom2 = groupList(jb)
865	>
866	>	if (skipThisPair(atom1, atom2)) cycle inner
867	>
868	>	call get_interatomic_vector(q(:,atom1), &
869	>	q(:,atom2), d_atm, ratmsq)
870	>
871	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
872	>	do_pot, do_stress, &
873	>	u_l, A, f, t, pot, vpair)
874	>
875	>	vab = vab + vpair
876	>
877	>	enddo inner
878	>	enddo
879	>
880	>	if (in_switching_region) then
881	>	swderiv = vab*dswdr/rgrp
882	>	do ia=groupStart(i), groupStart(i+1)-1
883	>	atom1=groupList(ia)
884	>	mf = mfact(atom1)
885	>	f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
886	>	f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
887	>	f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
888	>	enddo
889	>
890	>	do jb=groupStart(j), groupStart(j+1)-1
891	>	atom2=groupList(jb)
892	>	mf = mfact(atom2)
893	>	f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
894	>	f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
895	>	f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
896	>	enddo
897	>	endif
898	>
899	>	end if
900	>	enddo
901		enddo
902	+	point(nGroup) = nlist + 1
903
904	<	point(natoms) = nlist + 1
904	>	else  !! (of update_check)
905
724	–	else !! (update)
725	–
906		! use the list to find the neighbors
907	<	do i = 1, natoms-1
907	>	do i = 1, nGroup-1
908		JBEG = POINT(i)
909		JEND = POINT(i+1) - 1
910	<	! check thiat molecule i has neighbors
910	>	! check that group i has neighbors
911		if (jbeg .le. jend) then
912
913		do jnab = jbeg, jend
914		j = list(jnab)
915
916	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
917	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
738	<	u_l, A, f, t, pot)
916	>	call get_interatomic_vector(q_group(:,i), &
917	>	q_group(:,j), d_grp, rgrpsq)
918
919	+	vab = 0.0d0
920	+	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
921	+	in_switching_region)
922	+
923	+	do ia = groupStart(i), groupStart(i+1)-1
924	+	atom1 = groupList(ia)
925	+
926	+	do jb = groupStart(j), groupStart(j+1)-1
927	+	atom2 = groupList(jb)
928	+
929	+	call get_interatomic_vector(q(:,atom1), &
930	+	q(:,atom2), d_atm, ratmsq)
931	+
932	+	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
933	+	do_pot, do_stress, &
934	+	u_l, A, f, t, pot, vpair)
935	+
936	+	vab = vab + vpair
937	+
938	+	enddo
939	+	enddo
940	+
941	+	if (in_switching_region) then
942	+	swderiv = vab*dswdr/rgrp
943	+
944	+	do ia=groupStart(i), groupStart(i+1)-1
945	+	atom1=groupList(ia)
946	+	mf = mfact(atom1)
947	+	f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
948	+	f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
949	+	f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
950	+	enddo
951	+
952	+	do jb=groupStart(j), groupStart(j+1)-1
953	+	atom2=groupList(jb)
954	+	mf = mfact(atom2)
955	+	f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
956	+	f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
957	+	f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
958	+	enddo
959	+	endif
960		enddo
961		endif
962		enddo
963		endif
964	<
964	>
965		#endif
966
967		! phew, done with main loop.
968	<
969	<	!! Do timing
968	>
969	>	!! Do timing
970		#ifdef PROFILE
971		call cpu_time(forceTimeFinal)
972		forceTime = forceTime + forceTimeFinal - forceTimeInitial
973	<	#endif
974	<
755	<
973	>	#endif
974	>
975		#ifdef IS_MPI
976		!!distribute forces
977	<
977	>
978		f_temp = 0.0_dp
979		call scatter(f_Row,f_temp,plan_row3d)
980		do i = 1,nlocal
981		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
982		end do
983	<
983	>
984		f_temp = 0.0_dp
985		call scatter(f_Col,f_temp,plan_col3d)
986		do i = 1,nlocal
#	Line 785 | Line 1004 | contains
1004		if (do_pot) then
1005		! scatter/gather pot_row into the members of my column
1006		call scatter(pot_Row, pot_Temp, plan_row)
1007	<
1007	>
1008		! scatter/gather pot_local into all other procs
1009		! add resultant to get total pot
1010		do i = 1, nlocal
#	Line 793 | Line 1012 | contains
1012		enddo
1013
1014		pot_Temp = 0.0_DP
1015	<
1015	>
1016		call scatter(pot_Col, pot_Temp, plan_col)
1017		do i = 1, nlocal
1018		pot_local = pot_local + pot_Temp(i)
1019		enddo
1020	<
1021	<	endif
1020	>
1021	>	endif
1022		#endif
1023	<
1023	>
1024		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1025
1026		if (FF_uses_RF .and. SIM_uses_RF) then
1027	<
1027	>
1028		#ifdef IS_MPI
1029		call scatter(rf_Row,rf,plan_row3d)
1030		call scatter(rf_Col,rf_Temp,plan_col3d)
#	Line 815 | Line 1034 | contains
1034		#endif
1035
1036		do i = 1, nLocal
1037	<
1037	>
1038		rfpot = 0.0_DP
1039		#ifdef IS_MPI
1040		me_i = atid_row(i)
1041		#else
1042		me_i = atid(i)
1043		#endif
1044	<
1044	>
1045		if (PropertyMap(me_i)%is_DP) then
1046	<
1046	>
1047		mu_i = PropertyMap(me_i)%dipole_moment
1048	<
1048	>
1049		!! The reaction field needs to include a self contribution
1050		!! to the field:
1051		call accumulate_self_rf(i, mu_i, u_l)
#	Line 843 | Line 1062 | contains
1062		enddo
1063		endif
1064		endif
1065	<
1066	<
1065	>
1066	>
1067		#ifdef IS_MPI
1068	<
1068	>
1069		if (do_pot) then
1070		pot = pot + pot_local
1071		!! we assume the c code will do the allreduce to get the total potential
1072		!! we could do it right here if we needed to...
1073		endif
1074	<
1074	>
1075		if (do_stress) then
1076	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1076	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1077		mpi_comm_world,mpi_err)
1078		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1079		mpi_comm_world,mpi_err)
1080		endif
1081	<
1081	>
1082		#else
1083	<
1083	>
1084		if (do_stress) then
1085		tau = tau_Temp
1086		virial = virial_Temp
#	Line 870 | Line 1089 | contains
1089		#endif
1090
1091
873	–
1092		end subroutine do_force_loop
1093
1094	<	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
1094	>
1095	>	subroutine do_pair(i, j, rijsq, d, sw, do_pot, do_stress, &
1096	>	u_l, A, f, t, pot, vpair)
1097
1098	<	real( kind = dp ) :: pot
1098	>	real( kind = dp ) :: pot, vpair, sw
1099	>	real( kind = dp ), dimension(nLocal)   :: mfact
1100		real( kind = dp ), dimension(3,nLocal) :: u_l
1101	<	real (kind=dp), dimension(9,nLocal) :: A
1102	<	real (kind=dp), dimension(3,nLocal) :: f
1103	<	real (kind=dp), dimension(3,nLocal) :: t
1101	>	real( kind = dp ), dimension(9,nLocal) :: A
1102	>	real( kind = dp ), dimension(3,nLocal) :: f
1103	>	real( kind = dp ), dimension(3,nLocal) :: t
1104
1105		logical, intent(inout) :: do_pot, do_stress
1106		integer, intent(in) :: i, j
1107	<	real ( kind = dp ), intent(inout)    :: rijsq
1107	>	real ( kind = dp ), intent(inout) :: rijsq
1108		real ( kind = dp )                :: r
1109		real ( kind = dp ), intent(inout) :: d(3)
1110		integer :: me_i, me_j
#	Line 904 | Line 1125 | contains
1125		if (FF_uses_LJ .and. SIM_uses_LJ) then
1126
1127		if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
1128	<	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1128	>	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1129	>	do_stress)
1130		endif
1131
1132		endif
#	Line 912 | Line 1134 | contains
1134		if (FF_uses_charges .and. SIM_uses_charges) then
1135
1136		if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
1137	<	call do_charge_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1137	>	call do_charge_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1138	>	do_stress)
1139		endif
1140
1141		endif
#	Line 920 | Line 1143 | contains
1143		if (FF_uses_dipoles .and. SIM_uses_dipoles) then
1144
1145		if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
1146	<	call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
1146	>	call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
1147		do_pot, do_stress)
1148		if (FF_uses_RF .and. SIM_uses_RF) then
1149		call accumulate_rf(i, j, r, u_l)
#	Line 933 | Line 1156 | contains
1156		if (FF_uses_Sticky .and. SIM_uses_sticky) then
1157
1158		if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1159	<	call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
1159	>	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, pot, A, f, t, &
1160		do_pot, do_stress)
1161		endif
1162
#	Line 943 | Line 1166 | contains
1166		if (FF_uses_GB .and. SIM_uses_GB) then
1167
1168		if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1169	<	call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
1169	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
1170		do_pot, do_stress)
1171		endif
1172
#	Line 952 | Line 1175 | contains
1175		if (FF_uses_EAM .and. SIM_uses_EAM) then
1176
1177		if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1178	<	call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1178	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, pot, f, &
1179	>	do_pot, do_stress)
1180		endif
1181
1182		endif
1183	<
1183	>
1184		end subroutine do_pair
1185
1186	<
1187	<
964	<	subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
1186	>	subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, &
1187	>	do_pot, do_stress, u_l, A, f, t, pot)
1188		real( kind = dp ) :: pot
1189		real( kind = dp ), dimension(3,nLocal) :: u_l
1190		real (kind=dp), dimension(9,nLocal) :: A
#	Line 970 | Line 1193 | contains
1193
1194		logical, intent(inout) :: do_pot, do_stress
1195		integer, intent(in) :: i, j
1196	<	real ( kind = dp ), intent(inout)    :: rijsq
1197	<	real ( kind = dp )                :: r
1198	<	real ( kind = dp ), intent(inout) :: d(3)
1196	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1197	>	real ( kind = dp )                :: r, rc
1198	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1199
1200		logical :: is_EAM_i, is_EAM_j
1201
1202		integer :: me_i, me_j
1203
1204	<	r = sqrt(rijsq)
1204	>
1205	>	r = sqrt(rijsq)
1206	>	if (SIM_uses_molecular_cutoffs) then
1207	>	rc = sqrt(rcijsq)
1208	>	else
1209	>	rc = r
1210	>	endif
1211
1212
1213		#ifdef IS_MPI
#	Line 1008 | Line 1237 | contains
1237
1238
1239
1240	<	subroutine do_preforce(nlocal,pot)
1241	<	integer :: nlocal
1242	<	real( kind = dp ) :: pot
1243	<
1244	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1245	<	call calc_EAM_preforce_Frho(nlocal,pot)
1246	<	endif
1247	<
1248	<
1249	<	end subroutine do_preforce
1250	<
1251	<
1252	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1253	<
1254	<	real (kind = dp), dimension(3) :: q_i
1255	<	real (kind = dp), dimension(3) :: q_j
1256	<	real ( kind = dp ), intent(out) :: r_sq
1257	<	real( kind = dp ) :: d(3), scaled(3)
1258	<	integer i
1259	<
1260	<	d(1:3) = q_j(1:3) - q_i(1:3)
1261	<
1262	<	! Wrap back into periodic box if necessary
1263	<	if ( SIM_uses_PBC ) then
1035	<
1036	<	if( .not.boxIsOrthorhombic ) then
1037	<	! calc the scaled coordinates.
1038	<
1039	<	scaled = matmul(HmatInv, d)
1040	<
1041	<	! wrap the scaled coordinates
1042	<
1043	<	scaled = scaled  - anint(scaled)
1044	<
1045	<
1046	<	! calc the wrapped real coordinates from the wrapped scaled
1047	<	! coordinates
1048	<
1049	<	d = matmul(Hmat,scaled)
1050	<
1051	<	else
1052	<	! calc the scaled coordinates.
1053	<
1054	<	do i = 1, 3
1055	<	scaled(i) = d(i) * HmatInv(i,i)
1056	<
1057	<	! wrap the scaled coordinates
1058	<
1059	<	scaled(i) = scaled(i) - anint(scaled(i))
1060	<
1061	<	! calc the wrapped real coordinates from the wrapped scaled
1062	<	! coordinates
1063	<
1064	<	d(i) = scaled(i)*Hmat(i,i)
1065	<	enddo
1066	<	endif
1067	<
1068	<	endif
1069	<
1070	<	r_sq = dot_product(d,d)
1071	<
1072	<	end subroutine get_interatomic_vector
1073	<
1074	<	subroutine zero_work_arrays()
1075	<
1076	<	#ifdef IS_MPI
1077	<
1078	<	q_Row = 0.0_dp
1079	<	q_Col = 0.0_dp
1080	<
1081	<	u_l_Row = 0.0_dp
1082	<	u_l_Col = 0.0_dp
1083	<
1084	<	A_Row = 0.0_dp
1085	<	A_Col = 0.0_dp
1086	<
1087	<	f_Row = 0.0_dp
1088	<	f_Col = 0.0_dp
1089	<	f_Temp = 0.0_dp
1240	>	subroutine do_preforce(nlocal,pot)
1241	>	integer :: nlocal
1242	>	real( kind = dp ) :: pot
1243	>
1244	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1245	>	call calc_EAM_preforce_Frho(nlocal,pot)
1246	>	endif
1247	>
1248	>
1249	>	end subroutine do_preforce
1250	>
1251	>
1252	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1253	>
1254	>	real (kind = dp), dimension(3) :: q_i
1255	>	real (kind = dp), dimension(3) :: q_j
1256	>	real ( kind = dp ), intent(out) :: r_sq
1257	>	real( kind = dp ) :: d(3), scaled(3)
1258	>	integer i
1259	>
1260	>	d(1:3) = q_j(1:3) - q_i(1:3)
1261	>
1262	>	! Wrap back into periodic box if necessary
1263	>	if ( SIM_uses_PBC ) then
1264
1265	<	t_Row = 0.0_dp
1266	<	t_Col = 0.0_dp
1267	<	t_Temp = 0.0_dp
1265	>	if( .not.boxIsOrthorhombic ) then
1266	>	! calc the scaled coordinates.
1267	>
1268	>	scaled = matmul(HmatInv, d)
1269	>
1270	>	! wrap the scaled coordinates
1271	>
1272	>	scaled = scaled  - anint(scaled)
1273	>
1274	>
1275	>	! calc the wrapped real coordinates from the wrapped scaled
1276	>	! coordinates
1277	>
1278	>	d = matmul(Hmat,scaled)
1279	>
1280	>	else
1281	>	! calc the scaled coordinates.
1282	>
1283	>	do i = 1, 3
1284	>	scaled(i) = d(i) * HmatInv(i,i)
1285	>
1286	>	! wrap the scaled coordinates
1287	>
1288	>	scaled(i) = scaled(i) - anint(scaled(i))
1289	>
1290	>	! calc the wrapped real coordinates from the wrapped scaled
1291	>	! coordinates
1292	>
1293	>	d(i) = scaled(i)*Hmat(i,i)
1294	>	enddo
1295	>	endif
1296	>
1297	>	endif
1298	>
1299	>	r_sq = dot_product(d,d)
1300	>
1301	>	end subroutine get_interatomic_vector
1302
1303	<	pot_Row = 0.0_dp
1304	<	pot_Col = 0.0_dp
1305	<	pot_Temp = 0.0_dp
1306	<
1307	<	rf_Row = 0.0_dp
1308	<	rf_Col = 0.0_dp
1101	<	rf_Temp = 0.0_dp
1102	<
1103	<	#endif
1303	>	subroutine zero_work_arrays()
1304	>
1305	>	#ifdef IS_MPI
1306	>
1307	>	q_Row = 0.0_dp
1308	>	q_Col = 0.0_dp
1309
1310	+	q_group_Row = 0.0_dp
1311	+	q_group_Col = 0.0_dp
1312	+
1313	+	u_l_Row = 0.0_dp
1314	+	u_l_Col = 0.0_dp
1315	+
1316	+	A_Row = 0.0_dp
1317	+	A_Col = 0.0_dp
1318	+
1319	+	f_Row = 0.0_dp
1320	+	f_Col = 0.0_dp
1321	+	f_Temp = 0.0_dp
1322	+
1323	+	t_Row = 0.0_dp
1324	+	t_Col = 0.0_dp
1325	+	t_Temp = 0.0_dp
1326	+
1327	+	pot_Row = 0.0_dp
1328	+	pot_Col = 0.0_dp
1329	+	pot_Temp = 0.0_dp
1330	+
1331	+	rf_Row = 0.0_dp
1332	+	rf_Col = 0.0_dp
1333	+	rf_Temp = 0.0_dp
1334	+
1335	+	#endif
1336
1337	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1338	<	call clean_EAM()
1339	<	endif
1340	<
1341	<
1342	<
1343	<
1344	<
1345	<	rf = 0.0_dp
1346	<	tau_Temp = 0.0_dp
1347	<	virial_Temp = 0.0_dp
1348	<	end subroutine zero_work_arrays
1349	<
1350	<	function skipThisPair(atom1, atom2) result(skip_it)
1351	<	integer, intent(in) :: atom1
1352	<	integer, intent(in), optional :: atom2
1353	<	logical :: skip_it
1354	<	integer :: unique_id_1, unique_id_2
1355	<	integer :: me_i,me_j
1356	<	integer :: i
1357	<
1358	<	skip_it = .false.
1359	<
1360	<	!! there are a number of reasons to skip a pair or a particle
1361	<	!! mostly we do this to exclude atoms who are involved in short
1131	<	!! range interactions (bonds, bends, torsions), but we also need
1132	<	!! to exclude some overcounted interactions that result from
1133	<	!! the parallel decomposition
1134	<
1337	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1338	>	call clean_EAM()
1339	>	endif
1340	>
1341	>	rf = 0.0_dp
1342	>	tau_Temp = 0.0_dp
1343	>	virial_Temp = 0.0_dp
1344	>	end subroutine zero_work_arrays
1345	>
1346	>	function skipThisPair(atom1, atom2) result(skip_it)
1347	>	integer, intent(in) :: atom1
1348	>	integer, intent(in), optional :: atom2
1349	>	logical :: skip_it
1350	>	integer :: unique_id_1, unique_id_2
1351	>	integer :: me_i,me_j
1352	>	integer :: i
1353	>
1354	>	skip_it = .false.
1355	>
1356	>	!! there are a number of reasons to skip a pair or a particle
1357	>	!! mostly we do this to exclude atoms who are involved in short
1358	>	!! range interactions (bonds, bends, torsions), but we also need
1359	>	!! to exclude some overcounted interactions that result from
1360	>	!! the parallel decomposition
1361	>
1362		#ifdef IS_MPI
1363	<	!! in MPI, we have to look up the unique IDs for each atom
1364	<	unique_id_1 = tagRow(atom1)
1363	>	!! in MPI, we have to look up the unique IDs for each atom
1364	>	unique_id_1 = tagRow(atom1)
1365		#else
1366	<	!! in the normal loop, the atom numbers are unique
1367	<	unique_id_1 = atom1
1366	>	!! in the normal loop, the atom numbers are unique
1367	>	unique_id_1 = atom1
1368		#endif
1369	<
1370	<	!! We were called with only one atom, so just check the global exclude
1371	<	!! list for this atom
1372	<	if (.not. present(atom2)) then
1373	<	do i = 1, nExcludes_global
1374	<	if (excludesGlobal(i) == unique_id_1) then
1375	<	skip_it = .true.
1376	<	return
1377	<	end if
1378	<	end do
1379	<	return
1380	<	end if
1381	<
1369	>
1370	>	!! We were called with only one atom, so just check the global exclude
1371	>	!! list for this atom
1372	>	if (.not. present(atom2)) then
1373	>	do i = 1, nExcludes_global
1374	>	if (excludesGlobal(i) == unique_id_1) then
1375	>	skip_it = .true.
1376	>	return
1377	>	end if
1378	>	end do
1379	>	return
1380	>	end if
1381	>
1382		#ifdef IS_MPI
1383	<	unique_id_2 = tagColumn(atom2)
1383	>	unique_id_2 = tagColumn(atom2)
1384		#else
1385	<	unique_id_2 = atom2
1385	>	unique_id_2 = atom2
1386		#endif
1387	<
1387	>
1388		#ifdef IS_MPI
1389	<	!! this situation should only arise in MPI simulations
1390	<	if (unique_id_1 == unique_id_2) then
1391	<	skip_it = .true.
1392	<	return
1393	<	end if
1394	<
1395	<	!! this prevents us from doing the pair on multiple processors
1396	<	if (unique_id_1 < unique_id_2) then
1397	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1398	<	skip_it = .true.
1399	<	return
1400	<	endif
1401	<	else
1402	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1403	<	skip_it = .true.
1404	<	return
1405	<	endif
1406	<	endif
1389	>	!! this situation should only arise in MPI simulations
1390	>	if (unique_id_1 == unique_id_2) then
1391	>	skip_it = .true.
1392	>	return
1393	>	end if
1394	>
1395	>	!! this prevents us from doing the pair on multiple processors
1396	>	if (unique_id_1 < unique_id_2) then
1397	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1398	>	skip_it = .true.
1399	>	return
1400	>	endif
1401	>	else
1402	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1403	>	skip_it = .true.
1404	>	return
1405	>	endif
1406	>	endif
1407		#endif
1408	+
1409	+	!! the rest of these situations can happen in all simulations:
1410	+	do i = 1, nExcludes_global
1411	+	if ((excludesGlobal(i) == unique_id_1) .or. &
1412	+	(excludesGlobal(i) == unique_id_2)) then
1413	+	skip_it = .true.
1414	+	return
1415	+	endif
1416	+	enddo
1417	+
1418	+	do i = 1, nExcludes_local
1419	+	if (excludesLocal(1,i) == unique_id_1) then
1420	+	if (excludesLocal(2,i) == unique_id_2) then
1421	+	skip_it = .true.
1422	+	return
1423	+	endif
1424	+	else
1425	+	if (excludesLocal(1,i) == unique_id_2) then
1426	+	if (excludesLocal(2,i) == unique_id_1) then
1427	+	skip_it = .true.
1428	+	return
1429	+	endif
1430	+	endif
1431	+	endif
1432	+	end do
1433	+
1434	+	return
1435	+	end function skipThisPair
1436
1437	<	!! the rest of these situations can happen in all simulations:
1438	<	do i = 1, nExcludes_global
1439	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1440	<	(excludesGlobal(i) == unique_id_2)) then
1441	<	skip_it = .true.
1442	<	return
1443	<	endif
1444	<	enddo
1445	<
1446	<	do i = 1, nExcludes_local
1447	<	if (excludesLocal(1,i) == unique_id_1) then
1448	<	if (excludesLocal(2,i) == unique_id_2) then
1449	<	skip_it = .true.
1450	<	return
1451	<	endif
1452	<	else
1198	<	if (excludesLocal(1,i) == unique_id_2) then
1199	<	if (excludesLocal(2,i) == unique_id_1) then
1200	<	skip_it = .true.
1201	<	return
1202	<	endif
1203	<	endif
1204	<	endif
1205	<	end do
1206	<
1207	<	return
1208	<	end function skipThisPair
1209	<
1210	<	function FF_UsesDirectionalAtoms() result(doesit)
1211	<	logical :: doesit
1212	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1213	<	FF_uses_GB .or. FF_uses_RF
1214	<	end function FF_UsesDirectionalAtoms
1215	<
1216	<	function FF_RequiresPrepairCalc() result(doesit)
1217	<	logical :: doesit
1218	<	doesit = FF_uses_EAM
1219	<	end function FF_RequiresPrepairCalc
1220	<
1221	<	function FF_RequiresPostpairCalc() result(doesit)
1222	<	logical :: doesit
1223	<	doesit = FF_uses_RF
1224	<	end function FF_RequiresPostpairCalc
1225	<
1437	>	function FF_UsesDirectionalAtoms() result(doesit)
1438	>	logical :: doesit
1439	>	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1440	>	FF_uses_GB .or. FF_uses_RF
1441	>	end function FF_UsesDirectionalAtoms
1442	>
1443	>	function FF_RequiresPrepairCalc() result(doesit)
1444	>	logical :: doesit
1445	>	doesit = FF_uses_EAM
1446	>	end function FF_RequiresPrepairCalc
1447	>
1448	>	function FF_RequiresPostpairCalc() result(doesit)
1449	>	logical :: doesit
1450	>	doesit = FF_uses_RF
1451	>	end function FF_RequiresPostpairCalc
1452	>
1453		#ifdef PROFILE
1454	<	function getforcetime() result(totalforcetime)
1455	<	real(kind=dp) :: totalforcetime
1456	<	totalforcetime = forcetime
1457	<	end function getforcetime
1454	>	function getforcetime() result(totalforcetime)
1455	>	real(kind=dp) :: totalforcetime
1456	>	totalforcetime = forcetime
1457	>	end function getforcetime
1458		#endif
1459	<
1460	<	!! This cleans componets of force arrays belonging only to fortran
1461	<
1459	>
1460	>	!! This cleans componets of force arrays belonging only to fortran
1461	>
1462		end module do_Forces

Diff Legend

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