Skip to content

Commit ec4435c

Browse files
craabreucraabreu
committed
Code clean-up after RESPA implementation
1 parent b636969 commit ec4435c

16 files changed

+153
-365
lines changed

‎Makefile

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ OBJECTS = $(call obj,EmDeeCode neighbor_lists EmDeeData ArBee math structs model
7070

7171
COMPUTES = $(addprefix compute_,$(MODELTERMS))
7272
ENERGYCOMPUTES = $(addprefix energy_compute_,pair coul)
73-
VIRIALCOMPUTES = $(addprefix virial_compute_,pair coul unsplit)
73+
VIRIALCOMPUTES = $(addprefix virial_compute_,pair coul)
7474

7575
TESTS = $(patsubst %.f90,%,$(wildcard $(TSTDIR)/*.f90))
7676

@@ -187,9 +187,6 @@ $(SRCDIR)/virial_compute_pair.f90: $(call src,$(PAIRMODELS)) $(SRCDIR)/make_viri
187187
$(SRCDIR)/virial_compute_coul.f90: $(call src,$(COULMODELS)) $(SRCDIR)/make_virial_compute.sh
188188
bash $(SRCDIR)/make_virial_compute.sh coul $(COULMODELS) > $@
189189

190-
$(SRCDIR)/virial_compute_unsplit.f90: $(call src,$(COULMODELS)) $(SRCDIR)/make_unsplit_compute.sh
191-
bash $(SRCDIR)/make_unsplit_compute.sh $(COULMODELS) > $@
192-
193190
$(OBJDIR)/models.o: $(call obj,$(ALLMODELS) $(addprefix modelClass_,$(MODELTERMS))) \
194191
$(call obj,$(KSPACE) modelClass_kspace) $(SRCDIR)/make_models_module.sh
195192
bash $(SRCDIR)/make_models_module.sh $(ALLMODELS) $(KSPACE) > $(SRCDIR)/models.f90

‎src/EmDeeCode.f90

Lines changed: 15 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,6 @@
2020
! TODO: Include coul and kspace terms in subroutine update_forces.
2121
! TODO: Discount long-range electrostatic terms of bonds, angles, and dihedrals
2222
! TODO: Change name of exported Julia wrapper
23-
! TODO: Add option of automatic versus manual force calculations
2423
! TODO: Replace type(c_ptr) arguments by array arguments when possible
2524
! TODO: Create indexing for having sequential body particles and free particles in arrays
2625

@@ -33,9 +32,9 @@ module EmDeeCode
3332

3433
private
3534

36-
character(11), parameter :: VERSION = "27 Feb 2017"
35+
character(11), parameter :: VERSION = "28 Feb 2017"
3736

38-
type, bind(C) :: tOpts
37+
type, bind(C), public :: tOpts
3938
logical(lb) :: Translate ! Flag to activate/deactivate translations
4039
logical(lb) :: Rotate ! Flag to activate/deactivate rotations
4140
integer(ib) :: RotationMode ! Algorithm used for free rotation of rigid bodies
@@ -45,7 +44,7 @@ module EmDeeCode
4544
real(rb) :: Alpha_P ! Momentum-multiplying constant in momentum equations
4645
end type tOpts
4746

48-
type, bind(C) :: tEnergy
47+
type, bind(C), public :: tEnergy
4948
real(rb) :: Potential ! Total potential energy of the system
5049
real(rb) :: Dispersion ! Dispersion (vdW) part of the potential energy
5150
real(rb) :: Coulomb ! Electrostatic part of the potential energy
@@ -59,7 +58,7 @@ module EmDeeCode
5958
logical(lb) :: UpToDate ! Flag to attest whether energies have been computed
6059
end type tEnergy
6160

62-
type, bind(C) :: tEmDee
61+
type, bind(C), public :: tEmDee
6362
integer(ib) :: Builds ! Number of neighbor list builds
6463
real(rb) :: PairTime ! Time taken in force calculations
6564
real(rb) :: TotalTime ! Total time since initialization
@@ -219,7 +218,9 @@ subroutine EmDee_switch_model_layer( md, layer ) bind(C,name="EmDee_switch_model
219218

220219
call c_f_pointer( md%data, me )
221220
if (layer /= me%layer) then
222-
if ((layer < 1).or.(layer > me%nlayers)) call error( "switch_model_layer", "out of range" )
221+
if ((layer < 1).or.(layer > me%nlayers)) then
222+
call error( "model layer switch", "selected layer is out of range" )
223+
end if
223224
if (me%initialized) call update_forces( md, layer )
224225
me%layer = layer
225226
me%other_layer = [(i,i=1,layer-1), (i,i=layer+1,me%nlayers)]
@@ -624,7 +625,7 @@ subroutine EmDee_upload( md, option, address ) bind(C,name="EmDee_upload")
624625
me%initialized = allocated( me%R )
625626
if (me%initialized) call perform_initialization( me, md%DOF, md%RotDoF )
626627
end if
627-
if (me%initialized) call update_forces
628+
if (me%initialized) call compute_all_forces( md )
628629

629630
case ("coordinates")
630631
if (.not.allocated( me%R )) allocate( me%R(3,me%natoms) )
@@ -636,7 +637,7 @@ subroutine EmDee_upload( md, option, address ) bind(C,name="EmDee_upload")
636637
me%initialized = me%Lbox > zero
637638
if (me%initialized) call perform_initialization( me, md%DOF, md%RotDoF )
638639
end if
639-
if (me%initialized) call update_forces
640+
if (me%initialized) call compute_all_forces( md )
640641

641642
case ("momenta")
642643
if (.not.me%initialized) call error( "upload", "box and coordinates have not been defined" )
@@ -687,11 +688,6 @@ subroutine assign_charges( thread, Qext )
687688
me%charged(first:last) = abs(Qext(first:last)) > epsilon(one)
688689
end subroutine assign_charges
689690
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
690-
subroutine update_forces
691-
call compute_all_forces( md )
692-
if (me%respa_active) call compute_fast_forces( md )
693-
end subroutine update_forces
694-
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
695691
subroutine upload( thread, origin, destination )
696692
integer, intent(in) :: thread
697693
real(rb), intent(in) :: origin(3,me%natoms)
@@ -972,6 +968,8 @@ subroutine EmDee_advance( md, dt ) bind(C,name="EmDee_advance")
972968

973969
allocate( F_slow(3,me%natoms) )
974970

971+
if (me%fast_required) call compute_fast_forces( md )
972+
975973
!$omp parallel num_threads(me%nthreads)
976974
block
977975
integer :: thread, i
@@ -999,6 +997,7 @@ subroutine EmDee_advance( md, dt ) bind(C,name="EmDee_advance")
999997
me%Lbox = CR*me%Lbox
1000998
me%InvL = one/me%Lbox
1001999
me%invL2 = me%invL*me%invL
1000+
me%Lbox3 = me%Lbox
10021001
end if
10031002

10041003
call compute_fast_forces( md )
@@ -1045,6 +1044,7 @@ subroutine EmDee_advance( md, dt ) bind(C,name="EmDee_advance")
10451044
me%Lbox = CR*me%Lbox
10461045
me%InvL = one/me%Lbox
10471046
me%invL2 = me%invL*me%invL
1047+
me%Lbox3 = me%Lbox
10481048
end if
10491049

10501050
call compute_all_forces( md )
@@ -1071,8 +1071,6 @@ subroutine EmDee_advance( md, dt ) bind(C,name="EmDee_advance")
10711071
md%Energy%Kinetic = sum(md%Energy%TransPart) + md%Energy%Rotational
10721072
end if
10731073

1074-
if (changeBox) me%Lbox3 = me%Lbox
1075-
10761074
end subroutine EmDee_advance
10771075

10781076
!===================================================================================================
@@ -1131,6 +1129,8 @@ subroutine compute_all_forces( md )
11311129
me%Energy = sum(me%threadEnergy,2)
11321130
end if
11331131

1132+
me%fast_required = me%respa_active
1133+
11341134
md%Energy%UpToDate = compute
11351135
time = omp_get_wtime()
11361136
md%pairTime = md%pairTime + time
@@ -1170,26 +1170,6 @@ subroutine compute_fast_forces( md )
11701170

11711171
end subroutine compute_fast_forces
11721172

1173-
!===================================================================================================
1174-
1175-
! subroutine compute_fast_forces( me )
1176-
! type(tData), intent(inout) :: me
1177-
1178-
! real(rb) :: Rs(3,me%natoms), Fs(3,me%natoms,me%nthreads)
1179-
1180-
! Rs = me%invL*me%R
1181-
! !$omp parallel num_threads(me%nthreads)
1182-
! block
1183-
! integer :: thread
1184-
! thread = omp_get_thread_num() + 1
1185-
! call compute_short_range_forces( me, thread, Rs, Fs(:,:,thread) )
1186-
! end block
1187-
! !$omp end parallel
1188-
! me%F = sum(Fs,3)
1189-
!! if (me%nbodies /= 0) call rigid_body_forces( me )
1190-
1191-
! end subroutine compute_fast_forces
1192-
11931173
!===================================================================================================
11941174

11951175
subroutine update_forces( md, layer )
@@ -1227,58 +1207,4 @@ end subroutine update_forces
12271207

12281208
!===================================================================================================
12291209

1230-
! function total_virial( md ) result( virial ) bind(C,name="total_virial")
1231-
! type(tEmDee), intent(in) :: md
1232-
! real(rb) :: virial
1233-
1234-
! integer :: ib, jb, ik, jk, i, j, itype, jtype
1235-
! real(rb) :: Fij(3), Rij(3), r2, invR2, invR, Wij, QiQj, WCij, rFc, fshift
1236-
! type(tData), pointer :: me
1237-
1238-
! call c_f_pointer( md%data, me )
1239-
! fshift = one/me%RcSq
1240-
1241-
! virial = zero
1242-
! do ib = 1, me%nbodies-1
1243-
! do jb = ib+1, me%nbodies
1244-
! Fij = zero
1245-
! do ik = 1, me%body(ib)%NP
1246-
! i = me%body(ib)%index(ik)
1247-
! itype = me%atomType(i)
1248-
! do jk = 1, me%body(jb)%NP
1249-
! j = me%body(jb)%index(jk)
1250-
! Rij = me%R(:,i) - me%R(:,j)
1251-
! Rij = Rij - me%Lbox*anint(me%invL*Rij)
1252-
! r2 = sum(Rij*Rij)
1253-
! if (r2 < me%RcSq) then
1254-
! invR2 = one/r2
1255-
! invR = sqrt(invR2)
1256-
! jtype = me%atomType(j)
1257-
! associate( pair => me%pair(itype,jtype,me%layer) )
1258-
! select type ( model => pair%model )
1259-
! include "virial_compute_pair.f90"
1260-
! end select
1261-
! if (pair%model%shifted_force) then
1262-
! rFc = pair%model%fshift/invR
1263-
! Wij = Wij - rFc
1264-
! end if
1265-
! if (me%charged(i).and.me%charged(j).and.pair%coulomb) then
1266-
! QiQj = pair%kCoul*me%charge(i)*me%charge(j)
1267-
! WCij = QiQj*(invR - fshift/invR)
1268-
! Wij = Wij + WCij
1269-
! end if
1270-
! end associate
1271-
! Fij = Fij + Wij*invR2*Rij
1272-
! end if
1273-
! end do
1274-
! end do
1275-
! Rij = me%body(ib)%rcm - me%body(jb)%rcm
1276-
! Rij = Rij - me%Lbox*anint(me%invL*Rij)
1277-
! Virial = Virial + sum(Fij*Rij)
1278-
! end do
1279-
! end do
1280-
! end function total_virial
1281-
1282-
!===================================================================================================
1283-
12841210
end module EmDeeCode

‎src/EmDeeData.f90

Lines changed: 22 additions & 87 deletions
Original file line numberDiff line numberDiff line change
@@ -121,9 +121,9 @@ module EmDeeData
121121
real(rb) :: fshift
122122
integer(ib) :: Npair ! Number of core-neighbor sweepings per MD step
123123
integer(ib) :: Nbond ! Number of bonded computations per core sweeping
124-
124+
logical :: fast_required
125125
real(rb), allocatable :: F_fast(:,:)
126-
type(pairContainer), allocatable :: shortPair(:,:,:)
126+
type(pairContainer), allocatable :: closePair(:,:,:)
127127

128128
end type tData
129129

@@ -377,12 +377,12 @@ subroutine perform_initialization( me, DoF, RotDoF )
377377
end if
378378

379379
if (me%respa_active) then
380-
me%shortPair = me%pair
380+
me%closePair = me%pair
381381
allocate( me%F_fast(3,me%natoms) )
382382
do i = 1, me%ntypes
383383
do j = 1, me%ntypes
384384
do layer = 1, me%nlayers
385-
associate ( pair => me%shortPair(i,j,layer)%model )
385+
associate ( pair => me%closePair(i,j,layer)%model )
386386
pair%shifted_force = .true.
387387
call pair % shifting_setup( me%respaRc )
388388
end associate
@@ -576,29 +576,6 @@ end function cross
576576
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
577577
end subroutine compute_dihedrals
578578

579-
!===================================================================================================
580-
581-
subroutine compute_pair_forces( me, thread, Rs, F, Wpair, Wcoul )
582-
type(tData), intent(inout) :: me
583-
integer, intent(in) :: thread
584-
real(rb), intent(in) :: Rs(3,me%natoms)
585-
real(rb), intent(out) :: F(3,me%natoms), Wpair, Wcoul
586-
587-
real(rb) :: Rc2
588-
589-
Rc2 = me%RcSq*me%invL2
590-
# include "compute.f90"
591-
592-
contains
593-
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
594-
elemental function pbc( x )
595-
real(rb), intent(in) :: x
596-
real(rb) :: pbc
597-
pbc = x - anint(x)
598-
end function pbc
599-
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
600-
end subroutine compute_pair_forces
601-
602579
!===================================================================================================
603580

604581
subroutine compute_pairs( me, thread, compute, Rs, F, Epair, Ecoul, Wpair, Wcoul )
@@ -638,61 +615,12 @@ subroutine compute_short_range_forces( me, thread, Rs, F )
638615
real(rb), intent(in) :: Rs(3,me%natoms)
639616
real(rb), intent(out) :: F(3,me%natoms)
640617

641-
integer :: i, j, k, m, itype, jtype, firstAtom, lastAtom
642-
real(rb) :: Rc2, r2, invR2, invR, Wij, Qi, QiQj, WCij, rFc
643-
real(rb) :: Rij(3), Ri(3), Fi(3), Fij(3)
644-
logical :: icharged, ijcharged
645-
real(rb), allocatable :: Rvec(:,:)
618+
real(rb) :: Rc2
646619

647620
Rc2 = me%respaRcSq*me%invL2
648-
F = zero
649-
associate ( neighbor => me%neighbor(thread) )
650-
firstAtom = me%cellAtom%first(me%threadCell%first(thread))
651-
lastAtom = me%cellAtom%last(me%threadCell%last(thread))
652-
do k = firstAtom, lastAtom
653-
i = me%cellAtom%item(k)
654-
itype = me%atomType(i)
655-
Qi = me%charge(i)
656-
icharged = me%charged(i)
657-
Ri = Rs(:,i)
658-
Fi = zero
659-
associate ( partner => me%shortPair(:,itype,me%layer), &
660-
jlist => neighbor%item(neighbor%first(i):neighbor%middle(i)) )
661-
Rvec = Rs(:,jlist)
662-
forall (m=1:size(jlist)) Rvec(:,m) = pbc(Ri - Rvec(:,m))
663-
do m = 1, size(jlist)
664-
j = jlist(m)
665-
Rij = Rvec(:,m)
666-
r2 = sum(Rij*Rij)
667-
if (r2 < Rc2) then
668-
invR2 = me%invL2/r2
669-
invR = sqrt(invR2)
670-
jtype = me%atomType(j)
671-
ijcharged = icharged.and.me%charged(j)
672-
associate( pair => partner(jtype) )
673-
select type ( model => pair%model )
674-
include "virial_compute_pair.f90"
675-
end select
676-
if (pair%model%shifted_force) then
677-
rFc = pair%model%fshift/invR
678-
Wij = Wij - rFc
679-
end if
680-
if (ijcharged.and.pair%coulomb) then
681-
QiQj = pair%kCoul*Qi*me%charge(j)
682-
WCij = QiQj*(invR - me%fshift/invR)
683-
Wij = Wij + WCij
684-
end if
685-
end associate
686-
Fij = Wij*invR2*Rij
687-
Fi = Fi + Fij
688-
F(:,j) = F(:,j) - Fij
689-
end if
690-
end do
691-
end associate
692-
F(:,i) = F(:,i) + Fi
693-
end do
694-
end associate
695-
F = me%Lbox*F
621+
# define fast
622+
# include "compute.f90"
623+
# undef fast
696624

697625
contains
698626
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -913,21 +841,28 @@ function rigid_body_virial( me ) result( Virial )
913841
type(tData), intent(in) :: me
914842
real(rb) :: Virial
915843

916-
!$omp parallel num_threads(me%nthreads) reduction(+:Virial)
917-
Virial = partial_body_virial( omp_get_thread_num() + 1 )
844+
real(rb) :: W(me%nthreads)
845+
846+
!$omp parallel num_threads(me%nthreads)
847+
block
848+
integer :: thread
849+
thread = omp_get_thread_num() + 1
850+
call partial_body_virial( thread, W(thread) )
851+
end block
918852
!$omp end parallel
853+
Virial = -sum(W)
919854

920855
contains
921856
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
922-
function partial_body_virial( thread ) result (Virial )
857+
subroutine partial_body_virial( thread, W )
923858
integer, intent(in) :: thread
924-
real(rb) :: Virial
859+
real(rb), intent(out) :: W
925860
integer :: i
926-
Virial = zero
861+
W = zero
927862
do i = (thread - 1)*me%threadBodies + 1, min(thread*me%threadBodies, me%nbodies)
928-
Virial = Virial - sum(me%F(:,me%body(i)%index)*me%body(i)%delta)
863+
W = W + sum(me%F(:,me%body(i)%index)*me%body(i)%delta)
929864
end do
930-
end function partial_body_virial
865+
end subroutine partial_body_virial
931866
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
932867
end function rigid_body_virial
933868

0 commit comments

Comments
 (0)