4 |
|
|
5 |
|
!! @author Charles F. Vardeman II |
6 |
|
!! @author Matthew Meineke |
7 |
< |
!! @version $Id: do_Forces.F90,v 1.17 2003-03-13 00:33:18 chuckv Exp $, $Date: 2003-03-13 00:33:18 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $ |
7 |
> |
!! @version $Id: do_Forces.F90,v 1.19 2003-03-13 17:45:54 gezelter Exp $, $Date: 2003-03-13 17:45:54 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $ |
8 |
|
|
9 |
|
|
10 |
|
|
16 |
|
use lj |
17 |
|
use sticky_pair |
18 |
|
use dipole_dipole |
19 |
+ |
use reaction_field |
20 |
|
|
21 |
|
#ifdef IS_MPI |
22 |
|
use mpiSimulation |
32 |
|
logical, save :: FF_uses_GB |
33 |
|
logical, save :: FF_uses_EAM |
34 |
|
|
34 |
– |
|
35 |
|
public :: init_FF |
36 |
|
public :: do_force_loop |
37 |
|
|
38 |
|
contains |
39 |
|
|
40 |
< |
subroutine init_FF(thisStat) |
40 |
> |
subroutine init_FF(LJ_mix_policy, use_RF_c, thisStat) |
41 |
> |
logical(kind = 2), intent(in) :: use_RF_c |
42 |
> |
logical :: use_RF_f |
43 |
|
integer, intent(out) :: thisStat |
44 |
< |
integer :: my_status |
45 |
< |
character(len = 100) :: mix_Policy |
44 |
> |
integer :: my_status, nMatches |
45 |
> |
character(len = 100) :: LJ_mix_Policy |
46 |
> |
integer, pointer :: MatchList(:) |
47 |
> |
|
48 |
> |
!! Fortran's version of a cast: |
49 |
> |
use_RF_f = use_RF_c |
50 |
|
|
51 |
< |
! be a smarter subroutine. |
46 |
< |
mix_Policy = "FIXME" |
51 |
> |
!! assume things are copacetic, unless they aren't |
52 |
|
thisStat = 0 |
53 |
< |
call init_lj_FF(mix_Policy,my_status) |
53 |
> |
|
54 |
> |
!! init_FF is called *after* all of the atom types have been |
55 |
> |
!! defined in atype_module using the new_atype subroutine. |
56 |
> |
!! |
57 |
> |
!! this will scan through the known atypes and figure out what |
58 |
> |
!! interactions are used by the force field. |
59 |
> |
|
60 |
> |
FF_uses_LJ = .false. |
61 |
> |
FF_uses_sticky = .false. |
62 |
> |
FF_uses_dipoles = .false. |
63 |
> |
FF_uses_GB = .false. |
64 |
> |
FF_uses_EAM = .false. |
65 |
> |
|
66 |
> |
call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList) |
67 |
> |
deallocate(MatchList) |
68 |
> |
if (nMatches .gt. 0) FF_uses_LJ = .true. |
69 |
> |
|
70 |
> |
call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList) |
71 |
> |
deallocate(MatchList) |
72 |
> |
if (nMatches .gt. 0) FF_uses_dipoles = .true. |
73 |
> |
|
74 |
> |
call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, & |
75 |
> |
MatchList) |
76 |
> |
deallocate(MatchList) |
77 |
> |
if (nMatches .gt. 0) FF_uses_Sticky = .true. |
78 |
> |
|
79 |
> |
call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList) |
80 |
> |
deallocate(MatchList) |
81 |
> |
if (nMatches .gt. 0) FF_uses_GB = .true. |
82 |
> |
|
83 |
> |
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
84 |
> |
deallocate(MatchList) |
85 |
> |
if (nMatches .gt. 0) FF_uses_EAM = .true. |
86 |
> |
|
87 |
> |
!! check to make sure the use_RF setting makes sense |
88 |
> |
if (use_RF_f) then |
89 |
> |
if (FF_uses_dipoles) then |
90 |
> |
FF_uses_RF = .true. |
91 |
> |
call initialize_rf() |
92 |
> |
else |
93 |
> |
write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' |
94 |
> |
thisStat = -1 |
95 |
> |
return |
96 |
> |
endif |
97 |
> |
endif |
98 |
> |
|
99 |
> |
call init_lj_FF(LJ_mix_Policy, my_status) |
100 |
|
if (my_status /= 0) then |
101 |
|
thisStat = -1 |
102 |
|
return |
212 |
|
|
213 |
|
!! save current configuration, construct neighbor list, |
214 |
|
!! and calculate forces |
215 |
< |
call save_neighborList(q) |
215 |
> |
call saveNeighborList(q) |
216 |
|
|
217 |
|
neighborListSize = getNeighborListSize() |
218 |
|
nlist = 0 |
222 |
|
|
223 |
|
inner: do j = 1, ncol |
224 |
|
|
225 |
< |
if (checkExcludes(i,j)) cycle inner |
225 |
> |
if (skipThisPair(i,j)) cycle inner |
226 |
|
|
227 |
|
call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) |
228 |
|
|
276 |
|
|
277 |
|
! save current configuration, contruct neighbor list, |
278 |
|
! and calculate forces |
279 |
< |
call save_neighborList(q) |
279 |
> |
call saveNeighborList(q) |
280 |
|
|
281 |
|
neighborListSize = getNeighborListSize() |
282 |
|
nlist = 0 |
286 |
|
|
287 |
|
inner: do j = i+1, natoms |
288 |
|
|
289 |
< |
if (checkExcludes(i,j)) cycle inner |
289 |
> |
if (skipThisPair(i,j)) cycle inner |
290 |
|
|
291 |
|
call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) |
292 |
|
|
442 |
|
|
443 |
|
end subroutine do_force_loop |
444 |
|
|
394 |
– |
|
395 |
– |
!! Calculate any pre-force loop components and update nlist if necessary. |
396 |
– |
subroutine do_preForce(updateNlist) |
397 |
– |
logical, intent(inout) :: updateNlist |
398 |
– |
|
399 |
– |
|
400 |
– |
|
401 |
– |
end subroutine do_preForce |
402 |
– |
|
403 |
– |
!! Calculate any post force loop components, i.e. reaction field, etc. |
404 |
– |
subroutine do_postForce() |
405 |
– |
|
406 |
– |
|
407 |
– |
|
408 |
– |
end subroutine do_postForce |
409 |
– |
|
445 |
|
subroutine do_pair(i, j, rijsq, d, do_pot, do_stress) |
446 |
|
|
447 |
|
real( kind = dp ) :: pot |
537 |
|
|
538 |
|
subroutine check_initialization(error) |
539 |
|
integer, intent(out) :: error |
540 |
< |
|
540 |
> |
|
541 |
|
error = 0 |
542 |
|
! Make sure we are properly initialized. |
543 |
< |
if (.not. do_Forces_initialized) then |
543 |
> |
if (.not. do_forces_initialized) then |
544 |
|
write(default_error,*) "ERROR: do_Forces has not been initialized!" |
545 |
|
error = -1 |
546 |
|
return |
547 |
|
endif |
548 |
+ |
|
549 |
|
#ifdef IS_MPI |
550 |
|
if (.not. isMPISimSet()) then |
551 |
|
write(default_error,*) "ERROR: mpiSimulation has not been initialized!" |
553 |
|
return |
554 |
|
endif |
555 |
|
#endif |
556 |
< |
|
556 |
> |
|
557 |
|
return |
558 |
|
end subroutine check_initialization |
559 |
|
|
583 |
|
pot_Col = 0.0_dp |
584 |
|
pot_Temp = 0.0_dp |
585 |
|
|
586 |
+ |
rf_Row = 0.0_dp |
587 |
+ |
rf_Col = 0.0_dp |
588 |
+ |
rf_Temp = 0.0_dp |
589 |
+ |
|
590 |
|
#endif |
591 |
|
|
592 |
+ |
rf = 0.0_dp |
593 |
|
tau_Temp = 0.0_dp |
594 |
|
virial_Temp = 0.0_dp |
595 |
|
|
596 |
|
end subroutine zero_work_arrays |
597 |
|
|
598 |
< |
|
599 |
< |
!! Function to properly build neighbor lists in MPI using newtons 3rd law. |
600 |
< |
!! We don't want 2 processors doing the same i j pair twice. |
601 |
< |
!! Also checks to see if i and j are the same particle. |
598 |
> |
function skipThisPair(atom1, atom2) result(skip_it) |
599 |
> |
|
600 |
> |
integer, intent(in) :: atom1 |
601 |
> |
integer, intent(in), optional :: atom2 |
602 |
> |
logical :: skip_it |
603 |
> |
integer :: unique_id_1, unique_id_2 |
604 |
> |
integer :: i |
605 |
|
|
606 |
< |
function checkExcludes(atom1,atom2) result(do_cycle) |
607 |
< |
!--------------- Arguments-------------------------- |
608 |
< |
! Index i |
609 |
< |
integer,intent(in) :: atom1 |
610 |
< |
! Index j |
611 |
< |
integer,intent(in), optional :: atom2 |
612 |
< |
! Result do_cycle |
569 |
< |
logical :: do_cycle |
570 |
< |
!--------------- Local variables-------------------- |
571 |
< |
integer :: tag_i |
572 |
< |
integer :: tag_j |
573 |
< |
integer :: i, j |
574 |
< |
!--------------- END DECLARATIONS------------------ |
575 |
< |
do_cycle = .false. |
606 |
> |
skip_it = .false. |
607 |
> |
|
608 |
> |
!! there are a number of reasons to skip a pair or a particle |
609 |
> |
!! mostly we do this to exclude atoms who are involved in short |
610 |
> |
!! range interactions (bonds, bends, torsions), but we also need |
611 |
> |
!! to exclude some overcounted interactions that result from |
612 |
> |
!! the parallel decomposition |
613 |
|
|
614 |
|
#ifdef IS_MPI |
615 |
< |
tag_i = tagRow(atom1) |
615 |
> |
!! in MPI, we have to look up the unique IDs for each atom |
616 |
> |
unique_id_1 = tagRow(atom1) |
617 |
|
#else |
618 |
< |
tag_i = tag(atom1) |
618 |
> |
!! in the normal loop, the atom numbers are unique |
619 |
> |
unique_id_1 = atom1 |
620 |
|
#endif |
621 |
|
|
622 |
< |
!! Check global excludes first |
622 |
> |
!! We were called with only one atom, so just check the global exclude |
623 |
> |
!! list for this atom |
624 |
|
if (.not. present(atom2)) then |
625 |
|
do i = 1, nExcludes_global |
626 |
< |
if (excludeGlobal(i) == tag_i) then |
627 |
< |
do_cycle = .true. |
626 |
> |
if (excludesGlobal(i) == unique_id_1) then |
627 |
> |
skip_it = .true. |
628 |
|
return |
629 |
|
end if |
630 |
|
end do |
631 |
< |
return !! return after checking globals |
631 |
> |
return |
632 |
|
end if |
593 |
– |
|
594 |
– |
!! we return if atom2 not present here. |
595 |
– |
tag_j = tagColumn(atom2) |
633 |
|
|
634 |
< |
if (tag_i == tag_j) then |
635 |
< |
do_cycle = .true. |
634 |
> |
#ifdef IS_MPI |
635 |
> |
unique_id_2 = tagColumn(atom2) |
636 |
> |
#else |
637 |
> |
unique_id_2 = atom2 |
638 |
> |
#endif |
639 |
> |
|
640 |
> |
#ifdef IS_MPI |
641 |
> |
!! this situation should only arise in MPI simulations |
642 |
> |
if (unique_id_1 == unique_id_2) then |
643 |
> |
skip_it = .true. |
644 |
|
return |
645 |
|
end if |
646 |
|
|
647 |
< |
if (tag_i < tag_j) then |
648 |
< |
if (mod(tag_i + tag_j,2) == 0) do_cycle = .true. |
647 |
> |
!! this prevents us from doing the pair on multiple processors |
648 |
> |
if (unique_id_1 < unique_id_2) then |
649 |
> |
if (mod(unique_id_1 + unique_id_2,2) == 0) skip_it = .true. |
650 |
|
return |
651 |
|
else |
652 |
< |
if (mod(tag_i + tag_j,2) == 1) do_cycle = .true. |
652 |
> |
if (mod(unique_id_1 + unique_id_2,2) == 1) skip_it = .true. |
653 |
|
endif |
654 |
< |
|
655 |
< |
do i = 1, nExcludes_local |
656 |
< |
if ((tag_i == excludesLocal(1,i)) .and. (excludesLocal(2,i) < 0)) then |
657 |
< |
do_cycle = .true. |
654 |
> |
#endif |
655 |
> |
|
656 |
> |
!! the rest of these situations can happen in all simulations: |
657 |
> |
do i = 1, nExcludes_global |
658 |
> |
if ((excludesGlobal(i) == unique_id_1) .or. & |
659 |
> |
(excludesGlobal(i) == unique_id_2)) then |
660 |
> |
skip_it = .true. |
661 |
|
return |
662 |
< |
end if |
662 |
> |
endif |
663 |
> |
enddo |
664 |
> |
|
665 |
> |
do i = 1, nExcludes_local |
666 |
> |
if (excludesLocal(1,i) == unique_id_1) then |
667 |
> |
if (excludesLocal(2,i) == unique_id_2) then |
668 |
> |
skip_it = .true. |
669 |
> |
return |
670 |
> |
endif |
671 |
> |
else |
672 |
> |
if (excludesLocal(1,i) == unique_id_2) then |
673 |
> |
if (excludesLocal(2,i) == unique_id_1) then |
674 |
> |
skip_it = .true. |
675 |
> |
return |
676 |
> |
endif |
677 |
> |
endif |
678 |
> |
endif |
679 |
|
end do |
680 |
|
|
681 |
< |
|
682 |
< |
end function checkExcludes |
681 |
> |
return |
682 |
> |
end function skipThisPair |
683 |
|
|
684 |
|
function FF_UsesDirectionalAtoms() result(doesit) |
685 |
|
logical :: doesit |