Worked on the large memory use

.gitignore

@@ -10,7 +10,7 @@ bin/*
 # Ignore produced files
 **/outputs/*
 **/*.dat
-**/log.maniac
+**/log*maniac
 
 # Ignore backup and temporary files (optional but useful)
 *~

run.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 set -e  # Exit on error
 
-case="ZIF8-CH4O"
+case="WIDOM"
 
 base_energy="mc-topology/testcase-energy"
 base_adsorption="mc-topology/testcase-adsorption"
@@ -16,10 +16,12 @@ case "$case" in
 
   "WIDOM")
     folder="$base_widom/ZIF8-MET"
+    # reservoir="$base_reservoir/CH4O"
     ;;
 
   "CO2")
     folder="$base_reservoir/CO2"
+    reservoir="$base_reservoir/CO2"
     ;;
 
   "SLIT")

src/data_parser.f90

@@ -692,7 +692,7 @@ subroutine read_lammps_bonds(infile, data_file_name, box)
 
         ! If no bonds are expected, skip the whole routine
         if (box%num%bonds == 0) then
-            call InfoMessage("No bonds expected in data file: " // trim(data_file_name))
+            call info_message("No bonds expected in data file: " // trim(data_file_name))
             return
         end if
 
@@ -702,7 +702,7 @@ subroutine read_lammps_bonds(infile, data_file_name, box)
 
             if (ios < 0) then
                 write(formatted_msg, '(A, A)') "No bonds found in data file: ", trim(data_file_name)
-                call InfoMessage(trim(formatted_msg))
+                call info_message(trim(formatted_msg))
             end if
 
             if (ios > 0) then
@@ -790,7 +790,7 @@ subroutine read_lammps_angles(infile, data_file_name, box)
         ! If no angles are expected, skip the whole routine
         if (box%num%angles == 0) then
             write(formatted_msg, '(A, A)') "No angles expected in data file: ", trim(data_file_name)
-            call InfoMessage(trim(formatted_msg))
+            call info_message(trim(formatted_msg))
             return
         end if
 
@@ -800,7 +800,7 @@ subroutine read_lammps_angles(infile, data_file_name, box)
 
             if (ios < 0) then
                 write(formatted_msg, '(A, A)') "No angles found in data file: ", trim(data_file_name)
-                call InfoMessage(trim(formatted_msg))
+                call info_message(trim(formatted_msg))
             end if
 
             if (ios > 0) then
@@ -895,7 +895,7 @@ subroutine read_lammps_dihedrals(infile, data_file_name, box)
 
             if (ios < 0) then
                 write(formatted_msg, '(A, A)') "No dihedrals found in data file: ", trim(data_file_name)
-                call InfoMessage(trim(formatted_msg))
+                call info_message(trim(formatted_msg))
             end if
 
             if (ios > 0) then
@@ -993,7 +993,7 @@ subroutine read_lammps_impropers(infile, data_file_name, box)
             read(infile, '(A)', IOSTAT=ios) line
             if (ios < 0) then
                 write(formatted_msg, '(A, A)') "No impropers found in data file: ", trim(data_file_name)
-                call InfoMessage(trim(formatted_msg))
+                call info_message(trim(formatted_msg))
             end if
             if (ios > 0) then
                 write(formatted_msg, '(A, A)') "I/O error reading data file: ", trim(data_file_name)
@@ -1286,7 +1286,7 @@ subroutine detect_molecules(box)
         ! Global check: maximum allowed molecules
         if (detected_nb_max_molecule > NB_MAX_MOLECULE) then
             write(formatted_msg, '(A, I0)') "The number of molecules exceeds the maximum allowed = ", NB_MAX_MOLECULE
-            call InfoMessage(trim(formatted_msg))
+            call info_message(trim(formatted_msg))
             write(formatted_msg, '(A, I0)') "Number of molecules found = ", detected_nb_max_molecule
             call abort_run(trim(formatted_msg) // " Please increase 'NB_MAX_MOLECULE' in src/parameters.f90 and recompile.", 11)
         end if

src/ewald_energy.f90

@@ -82,6 +82,9 @@ subroutine compute_recip_energy_single_mol(res_type, mol_index, u_recipCoulomb_n
         logical :: creation_flag
         logical :: deletion_flag
 
+        ! Ensure that the routine deals with guest, not host
+        if (res%role(res_type) == TYPE_HOST) call abort_run("Inconsistence in fourier routine")
+
         creation_flag = present_or_false(is_creation)
         deletion_flag = present_or_false(is_deletion)
 
@@ -105,9 +108,9 @@ subroutine compute_recip_energy_single_mol(res_type, mol_index, u_recipCoulomb_n
             kz_idx = ewald%kvectors(idx)%kz
 
             ! Compute phase factors for the current residue
-            ewald%phase%new(1:natoms) = ewald%phase%factor(1, res_type, mol_index, 1:natoms, kx_idx) * &
-                ewald%phase%factor(2, res_type, mol_index, 1:natoms, ky_idx) * &
-                ewald%phase%factor(3, res_type, mol_index, 1:natoms, kz_idx)
+            ewald%phase%new(1:natoms) = ewald%phase%factor_guest(1, res_type, mol_index, 1:natoms, kx_idx) * &
+                ewald%phase%factor_guest(2, res_type, mol_index, 1:natoms, ky_idx) * &
+                ewald%phase%factor_guest(3, res_type, mol_index, 1:natoms, kz_idx)
 
             ewald%phase%old(1:natoms) = ewald%phase%factor_old(1, 1:natoms, kx_idx) * &
                 ewald%phase%factor_old(2, 1:natoms, ky_idx) * &
@@ -265,23 +268,69 @@ pure function compute_recip_amplitude(kx_idx, ky_idx, kz_idx) result(Ak)
 
         ! Loop over all residue types
         do res_type = 1, res%number
-            ! Loop over all molecules of this residue type
-            do mol_index = 1, primary%num%residues(res_type)
-                ! Loop over sites in molecule
-                do atom_index = 1, res%atom(res_type)
-
-                    ! Extract charge and phase product
-                    charges = primary%atoms%charges(res_type, atom_index)
-                    product = ewald%phase%factor(1, res_type, mol_index, atom_index, kx_idx) * &
-                            ewald%phase%factor(2, res_type, mol_index, atom_index, ky_idx) * &
-                            ewald%phase%factor(3, res_type, mol_index, atom_index, kz_idx)
-
-                    ! Accumulate contribution from this atom:
-                    ! charge * exp(i k · r) factor in x, y, z directions
-                    Ak = Ak + charges * product
-
-                end do
-            end do
+
+            ! For host
+            if (res%role(res_type) == TYPE_HOST) then
+
+                ! Loop over all molecules of this residue type
+                do mol_index = 1, primary%num%residues(res_type)
+                    ! Loop over sites in molecule
+                    do atom_index = 1, res%atom(res_type)
+
+                        ! Extract charge and phase product
+                        charges = primary%atoms%charges(res_type, atom_index)
+                        product = ewald%phase%factor_host(1, res_type, mol_index, atom_index, kx_idx) * &
+                                ewald%phase%factor_host(2, res_type, mol_index, atom_index, ky_idx) * &
+                                ewald%phase%factor_host(3, res_type, mol_index, atom_index, kz_idx)
+
+                        ! Accumulate contribution from this atom:
+                        ! charge * exp(i k · r) factor in x, y, z directions
+                        Ak = Ak + charges * product
+
+                    end do
+                end do 
+
+            ! For guest
+            else if (res%role(res_type) == TYPE_GUEST) then
+
+                ! Loop over all molecules of this residue type
+                do mol_index = 1, primary%num%residues(res_type)
+                    ! Loop over sites in molecule
+                    do atom_index = 1, res%atom(res_type)
+
+                        ! Extract charge and phase product
+                        charges = primary%atoms%charges(res_type, atom_index)
+                        product = ewald%phase%factor_guest(1, res_type, mol_index, atom_index, kx_idx) * &
+                                ewald%phase%factor_guest(2, res_type, mol_index, atom_index, ky_idx) * &
+                                ewald%phase%factor_guest(3, res_type, mol_index, atom_index, kz_idx)
+
+                        ! Accumulate contribution from this atom:
+                        ! charge * exp(i k · r) factor in x, y, z directions
+                        Ak = Ak + charges * product
+
+                    end do
+                end do 
+
+            end if
+
+            ! ! Loop over all molecules of this residue type
+            ! do mol_index = 1, primary%num%residues(res_type)
+            !     ! Loop over sites in molecule
+            !     do atom_index = 1, res%atom(res_type)
+
+            !         ! Extract charge and phase product
+            !         charges = primary%atoms%charges(res_type, atom_index)
+            !         product = ewald%phase%factor(1, res_type, mol_index, atom_index, kx_idx) * &
+            !                 ewald%phase%factor(2, res_type, mol_index, atom_index, ky_idx) * &
+            !                 ewald%phase%factor(3, res_type, mol_index, atom_index, kz_idx)
+
+            !         ! Accumulate contribution from this atom:
+            !         ! charge * exp(i k · r) factor in x, y, z directions
+            !         Ak = Ak + charges * product
+
+            !     end do
+            ! end do
+
         end do
 
     end function compute_recip_amplitude

src/ewald_phase.f90

@@ -14,11 +14,11 @@ module ewald_phase
     ! Saves the current Fourier-space phase factors and reciprocal amplitudes
     ! for a specific molecule or residue, enabling rollback after a rejected move.
     !--------------------------------------------------------------------
-    subroutine save_single_mol_fourier_terms(residue_type, molecule_index, symmetrize_x)
+    subroutine save_single_mol_fourier_terms(res_type, mol_index, symmetrize_x)
 
         ! Input arguments
-        integer, intent(in) :: residue_type       ! Residue type identifier
-        integer, intent(in) :: molecule_index     ! Molecule index to save
+        integer, intent(in) :: res_type       ! Residue type identifier
+        integer, intent(in) :: mol_index     ! Molecule index to save
         logical, intent(in), optional :: symmetrize_x ! Whether to save negative kx
 
         ! Local variables
@@ -32,24 +32,27 @@ subroutine save_single_mol_fourier_terms(residue_type, molecule_index, symmetriz
         do_sym = .false.
         if (present(symmetrize_x)) do_sym = symmetrize_x
 
+        ! Ensure that the routine deals with guest, not host
+        if (res%role(res_type) == TYPE_HOST) call abort_run("Inconsistence in fourier routine")
+
         !----------------------------------------------------------
         ! Save per-atom Fourier phase factors (IKX, IKY, IKZ)
         !----------------------------------------------------------
-        do atom_index = 1, res%atom(residue_type)
+        do atom_index = 1, res%atom(res_type)
 
             do dim = 1, 3
 
                 do k_idx = 0, ewald%param%kmax(dim)
 
                     ! Positive k
                     ewald%phase%factor_old(dim, atom_index, k_idx) = &
-                        ewald%phase%factor(dim, residue_type, molecule_index, atom_index, k_idx)
+                        ewald%phase%factor_guest(dim, res_type, mol_index, atom_index, k_idx)
 
                     ! Negative k (only if non-zero)
                     if (k_idx /= 0) then
                         if (dim /= 1 .or. do_sym) then
                             ewald%phase%factor_old(dim, atom_index, -k_idx) = &
-                                ewald%phase%factor(dim, residue_type, molecule_index, atom_index, -k_idx)
+                                ewald%phase%factor_guest(dim, res_type, mol_index, atom_index, -k_idx)
                         end if
                     end if
 
@@ -72,12 +75,12 @@ end subroutine save_single_mol_fourier_terms
     ! Restores previously saved Fourier-space phase factors and reciprocal
     ! amplitudes for a molecule or residue after a rejected move.
     !--------------------------------------------------------------------
-    subroutine restore_single_mol_fourier(residue_type, molecule_index, symmetrize_x)
+    subroutine restore_single_mol_fourier(res_type, mol_index, symmetrize_x)
 
         ! Input arguments
-        integer, intent(in) :: residue_type       ! Residue type index
-        integer, intent(in) :: molecule_index     ! Molecule index to restore
-        logical, intent(in), optional :: symmetrize_x ! Restore negative kx if true
+        integer, intent(in) :: res_type                 ! Residue type index
+        integer, intent(in) :: mol_index                ! Molecule index to restore
+        logical, intent(in), optional :: symmetrize_x   ! Restore negative kx if true
 
         ! Local variables
         integer :: atom_index_1                   ! Atom index
@@ -90,23 +93,24 @@ subroutine restore_single_mol_fourier(residue_type, molecule_index, symmetrize_x
         do_sym = .false.
         if (present(symmetrize_x)) do_sym = symmetrize_x
 
-        !----------------------------------------------------------
+        ! Ensure that the routine deals with guest, not host
+        if (res%role(res_type) == TYPE_HOST) call abort_run("Inconsistence in fourier routine")
+
         ! Restore per-atom Fourier phase factors (IKX, IKY, IKZ)
-        !----------------------------------------------------------
-        do atom_index_1 = 1, res%atom(residue_type)
+        do atom_index_1 = 1, res%atom(res_type)
 
             do dim = 1, 3
 
                 do k_idx = 0, ewald%param%kmax(dim)
 
                     ! Positive k
-                    ewald%phase%factor(dim, residue_type, molecule_index, atom_index_1, k_idx) = &
+                    ewald%phase%factor_guest(dim, res_type, mol_index, atom_index_1, k_idx) = &
                         ewald%phase%factor_old(dim, atom_index_1, k_idx)
 
                     ! Negative k
                     if (k_idx /= 0) then
                         if (dim /= 1 .or. do_sym) then
-                            ewald%phase%factor(dim, residue_type, molecule_index, atom_index_1, -k_idx) = &
+                            ewald%phase%factor_guest(dim, res_type, mol_index, atom_index_1, -k_idx) = &
                                 ewald%phase%factor_old(dim, atom_index_1, -k_idx)
                         end if
                     end if
@@ -128,10 +132,10 @@ end subroutine restore_single_mol_fourier
     ! Replaces the Fourier-space phase factors of one molecule (`index_1`)
     ! with those from another molecule (`index_2`) of the same residue type.
     !--------------------------------------------------------------------
-    subroutine replace_fourier_terms_single_mol(residue_type, index_1, index_2, symmetrize_x)
+    subroutine replace_fourier_terms_single_mol(res_type, index_1, index_2, symmetrize_x)
 
         ! Input arguments
-        integer, intent(in) :: residue_type     ! Residue type index
+        integer, intent(in) :: res_type     ! Residue type index
         integer, intent(in) :: index_1          ! Destination molecule index
         integer, intent(in) :: index_2          ! Source molecule index
         logical, intent(in), optional :: symmetrize_x   ! Copy negative kx if true
@@ -146,22 +150,25 @@ subroutine replace_fourier_terms_single_mol(residue_type, index_1, index_2, symm
         do_sym = .false.
         if (present(symmetrize_x)) do_sym = symmetrize_x
 
+        ! Ensure that the routine deals with guest, not host
+        if (res%role(res_type) == TYPE_HOST) call abort_run("Inconsistence in fourier routine")
+
         ! Restore IKX, IKY, IKZ from backups
-        do atom_index_1 = 1, res%atom(residue_type)
+        do atom_index_1 = 1, res%atom(res_type)
 
             do dim = 1, 3
 
                 do k_idx = 0, ewald%param%kmax(dim)
 
                     ! Always copy positive (and zero) k
-                    ewald%phase%factor(dim, residue_type, index_1, atom_index_1, k_idx) = &
-                        ewald%phase%factor(dim, residue_type, index_2, atom_index_1, k_idx)
+                    ewald%phase%factor_guest(dim, res_type, index_1, atom_index_1, k_idx) = &
+                        ewald%phase%factor_guest(dim, res_type, index_2, atom_index_1, k_idx)
 
                     ! Copy negative k only when allowed (when d=1, only if symmetrize_x is true, or when d=2,3)
                     if (k_idx /= 0) then
                         if (dim /= 1 .or. do_sym) then
-                            ewald%phase%factor(dim, residue_type, index_1, atom_index_1, -k_idx) = &
-                                ewald%phase%factor(dim, residue_type, index_2, atom_index_1, -k_idx)
+                            ewald%phase%factor_guest(dim, res_type, index_1, atom_index_1, -k_idx) = &
+                                ewald%phase%factor_guest(dim, res_type, index_2, atom_index_1, -k_idx)
                         end if
                     end if
 
@@ -178,17 +185,17 @@ end subroutine replace_fourier_terms_single_mol
     subroutine compute_all_fourier_terms()
 
         ! Local variables
-        integer :: residue_type
-        integer :: molecule_index
+        integer :: res_type
+        integer :: mol_index
 
         ! Loop over all residue types
-        do residue_type = 1, res%number
+        do res_type = 1, res%number
 
             ! Loop over all molecules of this residue type
-            do molecule_index = 1, primary%num%residues(residue_type)
+            do mol_index = 1, primary%num%residues(res_type)
 
                 ! Compute Fourier terms for a single molecule
-                call single_mol_fourier_terms(residue_type, molecule_index)
+                call single_mol_fourier_terms(res_type, mol_index)
 
             end do
         end do
@@ -230,11 +237,18 @@ subroutine single_mol_fourier_terms(res_type, mol_index)
             ! along each Cartesian direction. These factors will be used repeatedly
             ! in the reciprocal-space sum for the Ewald energy.
             do idim = 1, 3
+
                 call compute_phase_factor(ewald%phase%axis(:), local_phase(idim), ewald%param%kmax(idim))
-                ewald%phase%factor(idim, res_type, mol_index, atom_index_1, &
-                    -ewald%param%kmax(idim):ewald%param%kmax(idim)) = ewald%phase%axis(:)
-            end do
 
+                ! Differentiate between host and guest
+                if (res%role(res_type) == TYPE_HOST) then                    
+                    ewald%phase%factor_host(idim, res_type, mol_index, atom_index_1, &
+                        -ewald%param%kmax(idim):ewald%param%kmax(idim)) = ewald%phase%axis(:)
+                else if (res%role(res_type) == TYPE_GUEST) then
+                    ewald%phase%factor_guest(idim, res_type, mol_index, atom_index_1, &
+                        -ewald%param%kmax(idim):ewald%param%kmax(idim)) = ewald%phase%axis(:)
+                end if
+            end do
         end do
 
     end subroutine single_mol_fourier_terms