From 6aa00f3a1d10cc4c1239c9867ba2df2fe12f69c7 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sat, 29 Nov 2025 19:52:37 +0100 Subject: [PATCH 01/13] Worked on memory use --- src/ewald_energy.f90 | 92 +++++++++++++++++++++++++++++++--------- src/main.f90 | 8 ++-- src/prepare_utils.f90 | 26 +++++++++--- src/simulation_state.f90 | 2 + 4 files changed, 97 insertions(+), 31 deletions(-) diff --git a/src/ewald_energy.f90 b/src/ewald_energy.f90 index dafcf6d..e4c14b6 100644 --- a/src/ewald_energy.f90 +++ b/src/ewald_energy.f90 @@ -105,9 +105,15 @@ 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) + + ! #to remove + ! 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) + ! #to remove 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 +271,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 diff --git a/src/main.f90 b/src/main.f90 index 1b1acc2..d10da8a 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -26,10 +26,10 @@ program MANIAC call setup_simulation_parameters() ! Set up MC parameters, initial checks call compute_system_energy(primary) ! Compute initial total energy for the main box - ! Step 4 :Monte Carlo simulation - call monte_carlo_loop() ! Main MC loop + ! ! Step 4 :Monte Carlo simulation + ! call monte_carlo_loop() ! Main MC loop - ! Step 5 :Final reporting and cleanup - call energy_report(final=.true.) ! Print energy and statistics + ! ! Step 5 :Final reporting and cleanup + ! call energy_report(final=.true.) ! Print energy and statistics end program MANIAC diff --git a/src/prepare_utils.f90 b/src/prepare_utils.f90 index 06aaec4..92bec94 100644 --- a/src/prepare_utils.f90 +++ b/src/prepare_utils.f90 @@ -62,18 +62,28 @@ subroutine allocate_array() allocate(ewald%Ak_old(1:ewald%param%nkvec)) allocate(ewald%form_factor(ewald%param%nkvec)) ! Note, there is no need for such as large vector + write (*,*) "NB_MAX_MOLECULE", NB_MAX_MOLECULE + write (*,*) "nmax%active_residues", nmax%active_residues + write (*,*) "nmax%inactive_residues", nmax%inactive_residues + write (*,*) "nmax%atoms_per_residue", nmax%atoms_per_residue + write (*,*) "nmax%atoms_active_residue", nmax%atoms_active_residue + write (*,*) "nmax%atoms_inactive_residue", nmax%atoms_inactive_residue + ! Allocate complex arrays for wave vector components kmax_max = maxval(ewald%param%kmax) - allocate(ewald%phase%factor(3, res%number, 0:NB_MAX_MOLECULE, 1:nmax%atoms_per_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? - allocate(ewald%phase%factor_old(3, 1:nmax%atoms_per_residue, -kmax_max:kmax_max)) + allocate(ewald%phase%factor_host(3, res%number, 0:nmax%inactive_residues, 1:nmax%atoms_inactive_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? + allocate(ewald%phase%factor_guest(3, res%number, 0:nmax%active_residues, 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? + allocate(ewald%phase%factor_old(3, 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) + + ! #tofix - remove + allocate(ewald%phase%factor(3, res%number, 0:nmax%active_residues, 1:nmax%atoms_per_residue, -kmax_max:kmax_max)) + !allocate(ewald%phase%factor(3, res%number, 0:NB_MAX_MOLECULE, 1:nmax%atoms_per_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? bytes = size(ewald%phase%factor, kind=8) * storage_size(ewald%phase%factor) / 8 gigabytes = real(bytes, kind=8) / (1024.0d0**3) - - if (gigabytes > 1.0d0) then - write(msg,'(A,I0,A,F6.2,A)') "Large allocation detected: ", bytes, " bytes (", gigabytes, " GB)" - call warn_user(trim(msg)) - end if + write(msg,'(A,I0,A,F6.2,A)') "Large allocation detected: ", bytes, " bytes (", gigabytes, " GB)" + call warn_user(trim(msg)) + ! #tofix - remove ! Allocate temporary arrays once allocate(ewald%phase%axis(-kmax_max:kmax_max)) @@ -95,6 +105,8 @@ subroutine allocate_array() ewald%phase%new = zero ewald%phase%old = zero ewald%q_buffer = zero + ewald%phase%factor_guest = 0 + ewald%phase%factor_host = 0 ! Allocate widom if (proba%widom > 0) then diff --git a/src/simulation_state.f90 b/src/simulation_state.f90 index 01db258..ff6d83a 100644 --- a/src/simulation_state.f90 +++ b/src/simulation_state.f90 @@ -257,6 +257,8 @@ module simulation_state ! Allocatable array for phase calculation !--------------------------------------------------------------------------- type phase_ewald + complex(real64), allocatable :: factor_host(:,:,:,:,:) ! 5D array of precomputed complex exponentials e^(i k·r) for reciprocal-space summations + complex(real64), allocatable :: factor_guest(:,:,:,:,:) ! 5D array of precomputed complex exponentials e^(i k·r) for reciprocal-space summations complex(real64), allocatable :: factor(:,:,:,:,:) ! 5D array of precomputed complex exponentials e^(i k·r) for reciprocal-space summations complex(real64), allocatable :: factor_old(:,:,:) ! Previous configuration phase factors used for updating A(k) complex(real64), allocatable :: new(:) ! Temporary 1D array of phase factors for the new configuration From f9844b95d18ddce990db107005d7ef2fb2295a53 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sat, 29 Nov 2025 21:06:26 +0100 Subject: [PATCH 02/13] Reduced memory use --- src/data_parser.f90 | 14 ++++---- src/ewald_energy.f90 | 9 ++--- src/ewald_phase.f90 | 74 +++++++++++++++++++++++---------------- src/main.f90 | 8 ++--- src/output_utils.f90 | 4 +-- src/parameters_parser.f90 | 2 +- src/prepare_utils.f90 | 21 ----------- src/simulation_state.f90 | 1 - 8 files changed, 61 insertions(+), 72 deletions(-) diff --git a/src/data_parser.f90 b/src/data_parser.f90 index 7f288c1..9594164 100644 --- a/src/data_parser.f90 +++ b/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 diff --git a/src/ewald_energy.f90 b/src/ewald_energy.f90 index e4c14b6..0284660 100644 --- a/src/ewald_energy.f90 +++ b/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) @@ -109,12 +112,6 @@ subroutine compute_recip_energy_single_mol(res_type, mol_index, u_recipCoulomb_n 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) - ! #to remove - ! 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) - ! #to remove - ewald%phase%old(1:natoms) = ewald%phase%factor_old(1, 1:natoms, kx_idx) * & ewald%phase%factor_old(2, 1:natoms, ky_idx) * & ewald%phase%factor_old(3, 1:natoms, kz_idx) diff --git a/src/ewald_phase.f90 b/src/ewald_phase.f90 index 39542aa..49ff4e6 100644 --- a/src/ewald_phase.f90 +++ b/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,10 +32,13 @@ 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 @@ -43,13 +46,13 @@ subroutine save_single_mol_fourier_terms(residue_type, molecule_index, symmetriz ! 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 diff --git a/src/main.f90 b/src/main.f90 index d10da8a..1b1acc2 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -26,10 +26,10 @@ program MANIAC call setup_simulation_parameters() ! Set up MC parameters, initial checks call compute_system_energy(primary) ! Compute initial total energy for the main box - ! ! Step 4 :Monte Carlo simulation - ! call monte_carlo_loop() ! Main MC loop + ! Step 4 :Monte Carlo simulation + call monte_carlo_loop() ! Main MC loop - ! ! Step 5 :Final reporting and cleanup - ! call energy_report(final=.true.) ! Print energy and statistics + ! Step 5 :Final reporting and cleanup + call energy_report(final=.true.) ! Print energy and statistics end program MANIAC diff --git a/src/output_utils.f90 b/src/output_utils.f90 index abd6725..24f2e09 100644 --- a/src/output_utils.f90 +++ b/src/output_utils.f90 @@ -622,14 +622,14 @@ end subroutine warn_user !----------------------------------------------------------------------------- ! Print a standardized informational message !----------------------------------------------------------------------------- - subroutine InfoMessage(info_msg) + subroutine info_message(info_msg) character(len=*), intent(in) :: info_msg ! Information description call log_message("INFO: " // trim(info_msg)) call flush(output_unit) - end subroutine InfoMessage + end subroutine info_message !----------------------------------------------------------------------------- ! Validates that the given molecule index does not exceed the maximum allowed diff --git a/src/parameters_parser.f90 b/src/parameters_parser.f90 index cd3c29c..0ebdc8f 100644 --- a/src/parameters_parser.f90 +++ b/src/parameters_parser.f90 @@ -149,7 +149,7 @@ subroutine apply_lorentz_berthelot() ! Print the "rule enforcement" warning only once globally if (.not. any(warned)) then - call InfoMessage("Enforcing the Lorentz-Berthelot rule") + call info_message("Enforcing the Lorentz-Berthelot rule") call log_message("typei typej epsilon sigma") end if diff --git a/src/prepare_utils.f90 b/src/prepare_utils.f90 index 92bec94..43d8e40 100644 --- a/src/prepare_utils.f90 +++ b/src/prepare_utils.f90 @@ -49,9 +49,6 @@ subroutine allocate_array() ! Local variable integer :: kmax_max - integer(kind=8) :: bytes - real(8) :: gigabytes - character(len=256) :: msg allocate(saved%offset(3, nmax%atoms_per_residue)) allocate(saved%com(3)) @@ -62,29 +59,12 @@ subroutine allocate_array() allocate(ewald%Ak_old(1:ewald%param%nkvec)) allocate(ewald%form_factor(ewald%param%nkvec)) ! Note, there is no need for such as large vector - write (*,*) "NB_MAX_MOLECULE", NB_MAX_MOLECULE - write (*,*) "nmax%active_residues", nmax%active_residues - write (*,*) "nmax%inactive_residues", nmax%inactive_residues - write (*,*) "nmax%atoms_per_residue", nmax%atoms_per_residue - write (*,*) "nmax%atoms_active_residue", nmax%atoms_active_residue - write (*,*) "nmax%atoms_inactive_residue", nmax%atoms_inactive_residue - ! Allocate complex arrays for wave vector components kmax_max = maxval(ewald%param%kmax) allocate(ewald%phase%factor_host(3, res%number, 0:nmax%inactive_residues, 1:nmax%atoms_inactive_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? allocate(ewald%phase%factor_guest(3, res%number, 0:nmax%active_residues, 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? allocate(ewald%phase%factor_old(3, 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) - ! #tofix - remove - allocate(ewald%phase%factor(3, res%number, 0:nmax%active_residues, 1:nmax%atoms_per_residue, -kmax_max:kmax_max)) - !allocate(ewald%phase%factor(3, res%number, 0:NB_MAX_MOLECULE, 1:nmax%atoms_per_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? - - bytes = size(ewald%phase%factor, kind=8) * storage_size(ewald%phase%factor) / 8 - gigabytes = real(bytes, kind=8) / (1024.0d0**3) - write(msg,'(A,I0,A,F6.2,A)') "Large allocation detected: ", bytes, " bytes (", gigabytes, " GB)" - call warn_user(trim(msg)) - ! #tofix - remove - ! Allocate temporary arrays once allocate(ewald%phase%axis(-kmax_max:kmax_max)) allocate(ewald%phase%new(nmax%atoms_per_residue)) @@ -99,7 +79,6 @@ subroutine allocate_array() ewald%Ak = zero ewald%Ak_old = zero ewald%form_factor = zero - ewald%phase%factor = zero ewald%phase%factor_old = zero ewald%phase%axis = zero ewald%phase%new = zero diff --git a/src/simulation_state.f90 b/src/simulation_state.f90 index ff6d83a..a0dd9f0 100644 --- a/src/simulation_state.f90 +++ b/src/simulation_state.f90 @@ -259,7 +259,6 @@ module simulation_state type phase_ewald complex(real64), allocatable :: factor_host(:,:,:,:,:) ! 5D array of precomputed complex exponentials e^(i k·r) for reciprocal-space summations complex(real64), allocatable :: factor_guest(:,:,:,:,:) ! 5D array of precomputed complex exponentials e^(i k·r) for reciprocal-space summations - complex(real64), allocatable :: factor(:,:,:,:,:) ! 5D array of precomputed complex exponentials e^(i k·r) for reciprocal-space summations complex(real64), allocatable :: factor_old(:,:,:) ! Previous configuration phase factors used for updating A(k) complex(real64), allocatable :: new(:) ! Temporary 1D array of phase factors for the new configuration complex(real64), allocatable :: old(:) ! Temporary 1D array of phase factors for the old configuration From 55ddac91257b5df9e91d40a20297938b4cf6e823 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sat, 29 Nov 2025 21:54:31 +0100 Subject: [PATCH 03/13] Fixed wrong number of moelcules in memory --- run.sh | 3 ++- src/prescan_files.f90 | 14 +++++++------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/run.sh b/run.sh index afefc3e..3dc33f6 100755 --- a/run.sh +++ b/run.sh @@ -1,7 +1,7 @@ #!/bin/bash set -e # Exit on error -case="ZIF8-CH4O" +case="CO2" base_energy="mc-topology/testcase-energy" base_adsorption="mc-topology/testcase-adsorption" @@ -20,6 +20,7 @@ case "$case" in "CO2") folder="$base_reservoir/CO2" + reservoir="$base_reservoir/CO2" ;; "SLIT") diff --git a/src/prescan_files.f90 b/src/prescan_files.f90 index 9a5a156..d4afb28 100644 --- a/src/prescan_files.f90 +++ b/src/prescan_files.f90 @@ -38,9 +38,7 @@ subroutine prescan_inputs() ! Pre-scan input files call prescan_input_file(path%input) ! Detect residue definitions and sizes call prescan_topology(path%topology, is_reservoir = .false.) ! Count residues in primary topology - if (status%reservoir_provided) then - call prescan_topology(path%reservoir, is_reservoir = .true.) ! Count residues in reservoir topology - end if + if (status%reservoir_provided) call prescan_topology(path%reservoir, is_reservoir = .true.) ! Count residues in reservoir topology end subroutine prescan_inputs @@ -107,13 +105,16 @@ subroutine predetect_number_info(infile) ! Detect start/end of residue block select case (trim(line)) case ('begin_residue') + in_residue_block = .true. cycle + case ('end_residue') + in_residue_block = .false. res%number = res%number + 1 - cycle + end select if (.not. in_residue_block) cycle @@ -217,6 +218,7 @@ subroutine predetect_type(infile) allocate(res_infos(res%number)) do res_id = 1, res%number + res_infos(res_id)%nb_res = 0 res_infos(res_id)%nb_types = 0 res_infos(res_id)%nb_atoms = 0 if (allocated(res_infos(res_id)%types)) deallocate(res_infos(res_id)%types) @@ -390,8 +392,6 @@ subroutine prescan_topology(filename, is_reservoir) !---------------------------------------------------------- do res_id = 1, res%number - res_infos(res_id)%nb_res = 0 - if (res_infos(res_id)%nb_atoms <= 0) cycle ! Skip invalid residue ! Compute main/reservoir counts as real then round @@ -414,7 +414,7 @@ subroutine prescan_topology(filename, is_reservoir) ! Determine maximum active and inactive residues !---------------------------------------------------------- do res_id = 1, res%number - nb_residue = res_infos(res_id)%nb_res(1) + res_infos(res_id)%nb_res(2) + nb_residue = res_infos(res_id)%nb_res(1) + res_infos(res_id)%nb_res(2) if (res_infos(res_id)%is_active) then if (nb_residue > nmax%active_residues) then nmax%active_residues = nb_residue From 2a86a797e98bbb35ce73692c72560220d4ec8696 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sat, 29 Nov 2025 22:39:27 +0100 Subject: [PATCH 04/13] Keep looking for inconsistencies --- run.sh | 2 +- src/monte_carlo_utils.f90 | 20 +++++++------------- src/prescan_files.f90 | 6 +++++- src/widom.f90 | 5 ++--- 4 files changed, 15 insertions(+), 18 deletions(-) diff --git a/run.sh b/run.sh index 3dc33f6..a70769f 100755 --- a/run.sh +++ b/run.sh @@ -1,7 +1,7 @@ #!/bin/bash set -e # Exit on error -case="CO2" +case="ZIF8-CH4O" base_energy="mc-topology/testcase-energy" base_adsorption="mc-topology/testcase-adsorption" diff --git a/src/monte_carlo_utils.f90 b/src/monte_carlo_utils.f90 index d1447c8..aad8cb7 100644 --- a/src/monte_carlo_utils.f90 +++ b/src/monte_carlo_utils.f90 @@ -534,20 +534,14 @@ subroutine insert_and_orient_molecule(residue_type, molecule_index, rand_mol_ind ! Copy geometry from reservoir or rotate if no reservoir if (status%reservoir_provided) then - if (present(rand_mol_index)) then - ! Pick a random (and existing) molecule in the reservoir - call random_number(random_nmb) - rand_mol_index = int(random_nmb * reservoir%num%residues(residue_type)) + 1 + ! Pick a random (and existing) molecule in the reservoir + call random_number(random_nmb) - ! Copy site offsets from the chosen molecule - guest%offset(:, residue_type, molecule_index, 1:res%atom(residue_type)) = & - gas%offset(:, residue_type, rand_mol_index, 1:res%atom(residue_type)) - - else - - ! Error: reservoir exists but no output index provided - call abort_run("Rand mol index must be provided when reservoir exists.", 1) - end if + rand_mol_index = int(random_nmb * reservoir%num%residues(residue_type)) + 1 + + ! Copy site offsets from the chosen molecule + guest%offset(:, residue_type, molecule_index, 1:res%atom(residue_type)) = & + gas%offset(:, residue_type, rand_mol_index, 1:res%atom(residue_type)) else diff --git a/src/prescan_files.f90 b/src/prescan_files.f90 index d4afb28..f49654f 100644 --- a/src/prescan_files.f90 +++ b/src/prescan_files.f90 @@ -38,7 +38,11 @@ subroutine prescan_inputs() ! Pre-scan input files call prescan_input_file(path%input) ! Detect residue definitions and sizes call prescan_topology(path%topology, is_reservoir = .false.) ! Count residues in primary topology - if (status%reservoir_provided) call prescan_topology(path%reservoir, is_reservoir = .true.) ! Count residues in reservoir topology + if (status%reservoir_provided) then + call prescan_topology(path%reservoir, is_reservoir = .true.) ! Count residues in reservoir topology + else + nmax%active_residues = NB_MAX_MOLECULE + end if end subroutine prescan_inputs diff --git a/src/widom.f90 b/src/widom.f90 index 0ed2188..91a6c43 100644 --- a/src/widom.f90 +++ b/src/widom.f90 @@ -35,7 +35,6 @@ subroutine widom_trial(residue_type, molecule_index) ! Local variables real(real64) :: deltaU ! Energy difference - integer :: rand_mol_index ! Randomly selected molecule index from the reservoir for copying geometry call check_molecule_index(molecule_index) @@ -53,7 +52,7 @@ subroutine widom_trial(residue_type, molecule_index) call save_single_mol_fourier_terms(residue_type, molecule_index) ! Generate random insertion position within the simulation box - call insert_and_orient_molecule(residue_type, molecule_index, rand_mol_index) + call insert_and_orient_molecule(residue_type, molecule_index) ! Compute new energy call compute_new_energy(residue_type, molecule_index, is_creation = .true.) @@ -84,9 +83,9 @@ subroutine accumulate_widom_weight(residue_type, deltaU) if (weight > error) then ! Correspond to a success counter%widom(2) = counter%widom(2) + 1 + statistic%weight(residue_type) = statistic%weight(residue_type) + weight end if - statistic%weight(residue_type) = statistic%weight(residue_type) + weight statistic%sample(residue_type) = statistic%sample(residue_type) + 1 end subroutine accumulate_widom_weight From b4ccc0346f3e6156c71dad08e4c65745fdaef591 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sat, 29 Nov 2025 22:41:47 +0100 Subject: [PATCH 05/13] Fixed test --- tests/readers/prescan/test_prescan_files.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 506dbbd..3154ac9 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -18,7 +18,7 @@ program test_prescan_files path%parameters= trim(base) // "parameters.inc" ! Expectation in ZIF8-H2O system - expected_active = 3 + expected_active = 5000 expected_inactive = 1 ! Run prescan From 4805876ab44fb3d664367ee9bf70af18a4782624 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sat, 29 Nov 2025 23:15:32 +0100 Subject: [PATCH 06/13] Fixed error in insert_and_orient_molecule --- run.sh | 3 ++- src/monte_carlo_utils.f90 | 21 +++++++++++++------- src/output_utils.f90 | 2 +- src/widom.f90 | 1 + src/write_utils.f90 | 2 +- tests/readers/prescan/test_prescan_files.f90 | 6 +++--- 6 files changed, 22 insertions(+), 13 deletions(-) diff --git a/run.sh b/run.sh index a70769f..e22c3f3 100755 --- a/run.sh +++ b/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,6 +16,7 @@ case "$case" in "WIDOM") folder="$base_widom/ZIF8-MET" + # reservoir="$base_reservoir/CH4O" ;; "CO2") diff --git a/src/monte_carlo_utils.f90 b/src/monte_carlo_utils.f90 index aad8cb7..0acfda8 100644 --- a/src/monte_carlo_utils.f90 +++ b/src/monte_carlo_utils.f90 @@ -516,6 +516,7 @@ subroutine insert_and_orient_molecule(residue_type, molecule_index, rand_mol_ind logical :: full_rotation ! Flag indicating whether a full 360° random rotation should be applied real(real64) :: random_nmb ! Uniform random number in [0,1), used for random index selection real(real64) :: trial_pos(3) ! Random numbers for initial molecule position in the box + integer :: mol_index ! Decide whether to place a random COM if (present(place_random_com)) then @@ -537,11 +538,17 @@ subroutine insert_and_orient_molecule(residue_type, molecule_index, rand_mol_ind ! Pick a random (and existing) molecule in the reservoir call random_number(random_nmb) - rand_mol_index = int(random_nmb * reservoir%num%residues(residue_type)) + 1 - - ! Copy site offsets from the chosen molecule - guest%offset(:, residue_type, molecule_index, 1:res%atom(residue_type)) = & - gas%offset(:, residue_type, rand_mol_index, 1:res%atom(residue_type)) + if (present(place_random_com)) then + rand_mol_index = int(random_nmb * reservoir%num%residues(residue_type)) + 1 + ! Copy site offsets from the chosen molecule + guest%offset(:, residue_type, molecule_index, 1:res%atom(residue_type)) = & + gas%offset(:, residue_type, rand_mol_index, 1:res%atom(residue_type)) + else + mol_index = int(random_nmb * reservoir%num%residues(residue_type)) + 1 + ! Copy site offsets from the chosen molecule + guest%offset(:, residue_type, molecule_index, 1:res%atom(residue_type)) = & + gas%offset(:, residue_type, mol_index, 1:res%atom(residue_type)) + end if else @@ -580,7 +587,7 @@ end subroutine reject_creation_move ! using Widom particle insertion method, and optionally the ideal chemical ! potential (μ_ideal) for reporting purposes. !------------------------------------------------------------------------------ - subroutine CalculateExcessMu() + subroutine calculate_excess_mu() ! Local variables integer :: type_residue ! Residue type index @@ -618,7 +625,7 @@ subroutine CalculateExcessMu() end if end do - end subroutine CalculateExcessMu + end subroutine calculate_excess_mu !--------------------------------------------------------------------------- ! Physically removes a molecule from the simulation box by replacing it diff --git a/src/output_utils.f90 b/src/output_utils.f90 index 24f2e09..c9b5eb9 100644 --- a/src/output_utils.f90 +++ b/src/output_utils.f90 @@ -227,7 +227,7 @@ subroutine PrintStatus() end if if (proba%widom > 0) then - write(tmp,'("S(",I0,"/",I0,")")') counter%widom(2), counter%widom(1) + write(tmp,'("W(",I0,"/",I0,")")') counter%widom(2), counter%widom(1) move_msg = trim(move_msg)//" "//trim(tmp) end if diff --git a/src/widom.f90 b/src/widom.f90 index 91a6c43..7ff96dd 100644 --- a/src/widom.f90 +++ b/src/widom.f90 @@ -62,6 +62,7 @@ subroutine widom_trial(residue_type, molecule_index) ! Compute Boltzmann weight for Widom insertion deltaU = new%total - old%total + call accumulate_widom_weight(residue_type, deltaU) end subroutine widom_trial diff --git a/src/write_utils.f90 b/src/write_utils.f90 index c047bf3..7dd2071 100644 --- a/src/write_utils.f90 +++ b/src/write_utils.f90 @@ -204,7 +204,7 @@ subroutine write_dat_widom(nb_block, file_status) if (proba%widom > 0) then ! Compute excess and total chemical potentials - call CalculateExcessMu() + call calculate_excess_mu() ! Loop over residue types do type_residue = 1, res%number diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 3154ac9..75f10be 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -13,9 +13,9 @@ program test_prescan_files ! Setup paths base = "../../../mc-topology/testcase-energy/ZIF8-H2O/" - path%input = trim(base) // "input.maniac" - path%topology = trim(base) // "topology.data" - path%parameters= trim(base) // "parameters.inc" + path%input = trim(base) // "input.maniac" + path%topology = trim(base) // "topology.data" + path%parameters = trim(base) // "parameters.inc" ! Expectation in ZIF8-H2O system expected_active = 5000 From de811a0a1f3e46fca898b182065478069cec6b84 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 10:38:48 +0100 Subject: [PATCH 07/13] Improved reader prescan --- src/prescan_files.f90 | 38 ++++++++++++++--- tests/readers/prescan/run-test.sh | 5 +-- tests/readers/prescan/test_prescan_files.f90 | 44 ++++++++++++++++++-- 3 files changed, 74 insertions(+), 13 deletions(-) diff --git a/src/prescan_files.f90 b/src/prescan_files.f90 index f49654f..64fb44f 100644 --- a/src/prescan_files.f90 +++ b/src/prescan_files.f90 @@ -44,6 +44,20 @@ subroutine prescan_inputs() nmax%active_residues = NB_MAX_MOLECULE end if + + + + + + + + + + + + + + end subroutine prescan_inputs !--------------------------------------------------------------------------- @@ -336,7 +350,6 @@ subroutine prescan_topology(filename, is_reservoir) real(real64) :: q, x, y, z ! Charge and coordinate integer :: res_id ! Current residue index integer :: type_id ! Loop counter - integer :: nb_residue ! Local variable for number of residue logical :: found ! Logical indicator for atom counting integer, allocatable :: residue_atom_count(:) ! Temporary dynamic table for residue atom count @@ -391,9 +404,24 @@ subroutine prescan_topology(filename, is_reservoir) close(unit) - !---------------------------------------------------------- + call compute_max_residue(residue_atom_count, is_reservoir) + + end subroutine prescan_topology + + !----------------------------------------------------------------------------- + ! Compute max numbers of residues used for memory allocation. + !----------------------------------------------------------------------------- + subroutine compute_max_residue(residue_atom_count, is_reservoir) + + ! Input parameters + integer, intent(in) :: residue_atom_count(:) ! Temporary dynamic table for residue atom count + logical, intent(in) :: is_reservoir ! Flag: true if reservoir file + + ! Local variable + integer :: res_id ! Current residue index + integer :: nb_residue ! Local variable for number of residue + ! Compute number of residues - !---------------------------------------------------------- do res_id = 1, res%number if (res_infos(res_id)%nb_atoms <= 0) cycle ! Skip invalid residue @@ -414,9 +442,7 @@ subroutine prescan_topology(filename, is_reservoir) end do - !---------------------------------------------------------- ! Determine maximum active and inactive residues - !---------------------------------------------------------- do res_id = 1, res%number nb_residue = res_infos(res_id)%nb_res(1) + res_infos(res_id)%nb_res(2) if (res_infos(res_id)%is_active) then @@ -430,6 +456,6 @@ subroutine prescan_topology(filename, is_reservoir) end if end do - end subroutine prescan_topology + end subroutine compute_max_residue end module prescan_files diff --git a/tests/readers/prescan/run-test.sh b/tests/readers/prescan/run-test.sh index b6fc2d3..c958098 100755 --- a/tests/readers/prescan/run-test.sh +++ b/tests/readers/prescan/run-test.sh @@ -5,11 +5,11 @@ set -euo pipefail BUILD_DIR="../../../build" MODULE_DIR="../../../include" -# Where test executables should go +# Path to test executables BIN_DIR="./build" mkdir -p "$BIN_DIR" -# Collect Maniac object files, excluding main.o +# Collect Maniac object files (excluding main.o) shopt -s nullglob MANIAC_OBJS=() for obj in "$BUILD_DIR"/*.o; do @@ -18,7 +18,6 @@ for obj in "$BUILD_DIR"/*.o; do done shopt -u nullglob -# ----------------------------- # COMPILE & RUN TESTS for TEST_SRC in test_*.f90; do [[ ! -f "$TEST_SRC" ]] && { diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 75f10be..1718586 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -11,15 +11,51 @@ program test_prescan_files logical :: pass character(len=LENPATH) :: base + !--------------------------------------------------------------------------- + ! ZIF8-H2O system without reservoir + !--------------------------------------------------------------------------- + ! Setup paths base = "../../../mc-topology/testcase-energy/ZIF8-H2O/" path%input = trim(base) // "input.maniac" path%topology = trim(base) // "topology.data" path%parameters = trim(base) // "parameters.inc" - ! Expectation in ZIF8-H2O system - expected_active = 5000 - expected_inactive = 1 + expected_active = NB_MAX_MOLECULE ! In absence of reservoir, use NB_MAX_MOLECULE + expected_inactive = 1 ! One host is expected + + ! Run prescan + call prescan_inputs() + + ! Validate results + pass = (nmax%active_residues == expected_active) .and. & + (nmax%inactive_residues == expected_inactive) + + if (.not. pass) then + if (nmax%active_residues /= expected_active) then + print *, " active_residues expected=", expected_active, " got=", nmax%active_residues + end if + if (nmax%inactive_residues /= expected_inactive) then + print *, " inactive_residues expected=", expected_inactive, " got=", nmax%inactive_residues + end if + stop 1 + end if + + !--------------------------------------------------------------------------- + ! ZIF8-CH4O system with reservoir + !--------------------------------------------------------------------------- + + ! Setup paths + base = "../../../mc-topology/testcase-adsorption/ZIF8-CH4O/" + path%input = trim(base) // "input.maniac" + path%topology = trim(base) // "topology.data" + path%parameters = trim(base) // "parameters.inc" + base = "../../../mc-topology//molecule-reservoir/CH4O-H2O/" + path%reservoir = trim(base) // "topology.data" + status%reservoir_provided = .true. + + expected_active = 5+1500 ! 5 molecules in main, 1500 in reservoir + expected_inactive = 1 ! One host is expected ! Run prescan call prescan_inputs() @@ -29,7 +65,6 @@ program test_prescan_files (nmax%inactive_residues == expected_inactive) if (.not. pass) then - print *, "❌ Test FAILED!" if (nmax%active_residues /= expected_active) then print *, " active_residues expected=", expected_active, " got=", nmax%active_residues end if @@ -39,4 +74,5 @@ program test_prescan_files stop 1 end if + end program test_prescan_files From bc04558d953d7133e0933e031fb64cb1ad549f0c Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 11:05:13 +0100 Subject: [PATCH 08/13] Improved test --- src/prescan_files.f90 | 27 +----- tests/readers/prescan/test_prescan_files.f90 | 98 +++++++++----------- 2 files changed, 47 insertions(+), 78 deletions(-) diff --git a/src/prescan_files.f90 b/src/prescan_files.f90 index 64fb44f..b78d7a3 100644 --- a/src/prescan_files.f90 +++ b/src/prescan_files.f90 @@ -44,20 +44,6 @@ subroutine prescan_inputs() nmax%active_residues = NB_MAX_MOLECULE end if - - - - - - - - - - - - - - end subroutine prescan_inputs !--------------------------------------------------------------------------- @@ -356,15 +342,11 @@ subroutine prescan_topology(filename, is_reservoir) allocate(residue_atom_count(res%number)) residue_atom_count = 0 - ! --------------------------------------------------------- ! Open topology file - ! --------------------------------------------------------- open(newunit=unit, file=filename, status='old', action='read', iostat=ios) call check_IO_status(filename, ios) - ! --------------------------------------------------------- ! Scan file - ! --------------------------------------------------------- do read(unit,'(A)', iostat=ios) line if (ios /= 0) exit @@ -381,24 +363,21 @@ subroutine prescan_topology(filename, is_reservoir) cycle end if - ! Note, the atome type is expected to be in second position - ! Read atom-ID, molecule-ID, atom-type + ! Read line (the atom type is expected to be in second position) read(line,*,iostat=ios) atom_id, mol_id, atom_type, q, x, y, z if (ios /= 0) cycle - ! --------------------------------------------------------- ! Classify this atom according to which residue type it belongs to - ! --------------------------------------------------------- found = .false. do res_id = 1, res%number do type_id = 1, res_infos(res_id)%nb_types if (atom_type == res_infos(res_id)%types(type_id)) then residue_atom_count(res_id) = residue_atom_count(res_id) + 1 found = .true. - exit ! stop inner loop + exit end if end do - if (found) exit ! stop outer loop + if (found) exit end do end do diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 1718586..4c64527 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -1,78 +1,68 @@ program test_prescan_files use, intrinsic :: iso_fortran_env, only: real64 - use prescan_files - implicit none - character(len=256) :: topology_path, input_file, data_file, inc_file - integer :: expected_active, expected_inactive - logical :: pass - character(len=LENPATH) :: base + character(len=LENPATH) :: base_main, base_res !--------------------------------------------------------------------------- - ! ZIF8-H2O system without reservoir + ! Test 1 : ZIF8-H2O without reservoir !--------------------------------------------------------------------------- + base_main = "../../../mc-topology/testcase-energy/ZIF8-H2O/" - ! Setup paths - base = "../../../mc-topology/testcase-energy/ZIF8-H2O/" - path%input = trim(base) // "input.maniac" - path%topology = trim(base) // "topology.data" - path%parameters = trim(base) // "parameters.inc" + ! Expected : in absence of reservoir, use NB_MAX_MOLECULE + call run_test(base_main, "", .false., NB_MAX_MOLECULE, 1) - expected_active = NB_MAX_MOLECULE ! In absence of reservoir, use NB_MAX_MOLECULE - expected_inactive = 1 ! One host is expected - - ! Run prescan - call prescan_inputs() + !--------------------------------------------------------------------------- + ! Test 2 : ZIF8-CH4O system with reservoir + !--------------------------------------------------------------------------- + base_main = "../../../mc-topology/testcase-adsorption/ZIF8-CH4O/" + base_res = "../../../mc-topology/molecule-reservoir/CH4O-H2O/" - ! Validate results - pass = (nmax%active_residues == expected_active) .and. & - (nmax%inactive_residues == expected_inactive) + ! Expected : 5 molecules in main, 1500 in reservoir + call run_test(base_main, base_res, .true., 5+1500, 1) - if (.not. pass) then - if (nmax%active_residues /= expected_active) then - print *, " active_residues expected=", expected_active, " got=", nmax%active_residues - end if - if (nmax%inactive_residues /= expected_inactive) then - print *, " inactive_residues expected=", expected_inactive, " got=", nmax%inactive_residues - end if - stop 1 - end if +contains - !--------------------------------------------------------------------------- - ! ZIF8-CH4O system with reservoir - !--------------------------------------------------------------------------- + subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) - ! Setup paths - base = "../../../mc-topology/testcase-adsorption/ZIF8-CH4O/" - path%input = trim(base) // "input.maniac" - path%topology = trim(base) // "topology.data" - path%parameters = trim(base) // "parameters.inc" - base = "../../../mc-topology//molecule-reservoir/CH4O-H2O/" - path%reservoir = trim(base) // "topology.data" - status%reservoir_provided = .true. + ! Input parameters + character(len=*), intent(in) :: base_main + character(len=*), intent(in) :: base_res + logical, intent(in) :: reservoir + integer, intent(in) :: exp_active, exp_inactive - expected_active = 5+1500 ! 5 molecules in main, 1500 in reservoir - expected_inactive = 1 ! One host is expected + ! Local variable + logical :: pass - ! Run prescan - call prescan_inputs() + ! Set paths + path%input = trim(base_main) // "input.maniac" + path%topology = trim(base_main) // "topology.data" + path%parameters = trim(base_main) // "parameters.inc" - ! Validate results - pass = (nmax%active_residues == expected_active) .and. & - (nmax%inactive_residues == expected_inactive) + status%reservoir_provided = reservoir - if (.not. pass) then - if (nmax%active_residues /= expected_active) then - print *, " active_residues expected=", expected_active, " got=", nmax%active_residues + if (reservoir) then + path%reservoir = trim(base_res) // "topology.data" end if - if (nmax%inactive_residues /= expected_inactive) then - print *, " inactive_residues expected=", expected_inactive, " got=", nmax%inactive_residues + + ! Run prescan logic + call prescan_inputs() + + ! Check results + pass = (nmax%active_residues == exp_active) .and. & + (nmax%inactive_residues == exp_inactive) + + if (.not. pass) then + + print *, "Test FAILED" + print *, " expected active =", exp_active, " got=", nmax%active_residues + print *, " expected inactive =", exp_inactive, " got=", nmax%inactive_residues + stop 1 + end if - stop 1 - end if + end subroutine run_test end program test_prescan_files From 176e036b030bcd700f046b5bdb8a5a5994783b49 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 11:15:57 +0100 Subject: [PATCH 09/13] Improved test prescan --- mc-topology | 2 +- tests/readers/prescan/test_prescan_files.f90 | 26 ++++++++++++++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/mc-topology b/mc-topology index c0e6682..832e16b 160000 --- a/mc-topology +++ b/mc-topology @@ -1 +1 @@ -Subproject commit c0e66827a310c716478179fc224abdea730a281e +Subproject commit 832e16b24a80ad9211a254f3db123d43286ce4b6 diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 4c64527..5cd56e0 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -2,13 +2,16 @@ program test_prescan_files use, intrinsic :: iso_fortran_env, only: real64 use prescan_files + implicit none - character(len=LENPATH) :: base_main, base_res + character(len=LENPATH) :: base_main + character(len=LENPATH) :: base_res !--------------------------------------------------------------------------- ! Test 1 : ZIF8-H2O without reservoir !--------------------------------------------------------------------------- + base_main = "../../../mc-topology/testcase-energy/ZIF8-H2O/" ! Expected : in absence of reservoir, use NB_MAX_MOLECULE @@ -17,12 +20,31 @@ program test_prescan_files !--------------------------------------------------------------------------- ! Test 2 : ZIF8-CH4O system with reservoir !--------------------------------------------------------------------------- + base_main = "../../../mc-topology/testcase-adsorption/ZIF8-CH4O/" base_res = "../../../mc-topology/molecule-reservoir/CH4O-H2O/" ! Expected : 5 molecules in main, 1500 in reservoir call run_test(base_main, base_res, .true., 5+1500, 1) + !--------------------------------------------------------------------------- + ! Test 3 : Bulk CO2 system without reservoir + !--------------------------------------------------------------------------- + + base_main = "../../../mc-topology/molecule-reservoir/CO2/" + + ! Expected : in absence of reservoir, use NB_MAX_MOLECULE + call run_test(base_main, "", .false., NB_MAX_MOLECULE, 0) + + !--------------------------------------------------------------------------- + ! Test 4 : Bulk CO2 system with reservoir + !--------------------------------------------------------------------------- + + base_main = "../../../mc-topology/molecule-reservoir/CO2/" + + ! Expected : 50 molecules in main, 50 in reservoir + call run_test(base_main, base_main, .true., 100, 0) + contains subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) @@ -51,7 +73,7 @@ subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) call prescan_inputs() ! Check results - pass = (nmax%active_residues == exp_active) .and. & + pass = (nmax%active_residues == exp_active) .and. & (nmax%inactive_residues == exp_inactive) if (.not. pass) then From d54ec3c64e48fb3a35a3a60e156272705cdeea1c Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 11:20:51 +0100 Subject: [PATCH 10/13] Minor improvement of prescanned files --- tests/readers/prescan/test_prescan_files.f90 | 21 ++++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 5cd56e0..54f805b 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -13,9 +13,10 @@ program test_prescan_files !--------------------------------------------------------------------------- base_main = "../../../mc-topology/testcase-energy/ZIF8-H2O/" + base_res = "" ! Expected : in absence of reservoir, use NB_MAX_MOLECULE - call run_test(base_main, "", .false., NB_MAX_MOLECULE, 1) + call run_test(base_main, base_res, .false., NB_MAX_MOLECULE, 1) !--------------------------------------------------------------------------- ! Test 2 : ZIF8-CH4O system with reservoir @@ -32,9 +33,10 @@ program test_prescan_files !--------------------------------------------------------------------------- base_main = "../../../mc-topology/molecule-reservoir/CO2/" + base_res = "" ! Expected : in absence of reservoir, use NB_MAX_MOLECULE - call run_test(base_main, "", .false., NB_MAX_MOLECULE, 0) + call run_test(base_main, base_res, .false., NB_MAX_MOLECULE, 0) !--------------------------------------------------------------------------- ! Test 4 : Bulk CO2 system with reservoir @@ -50,10 +52,11 @@ program test_prescan_files subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) ! Input parameters - character(len=*), intent(in) :: base_main - character(len=*), intent(in) :: base_res - logical, intent(in) :: reservoir - integer, intent(in) :: exp_active, exp_inactive + character(len=*), intent(in) :: base_main ! Path for main file + character(len=*), intent(in) :: base_res ! Path for reservoir file (if reservoir is present) + logical, intent(in) :: reservoir ! Is reservoir present + integer, intent(in) :: exp_active ! Expected guest residue + integer, intent(in) :: exp_inactive ! Expected host residue ! Local variable logical :: pass @@ -65,9 +68,7 @@ subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) status%reservoir_provided = reservoir - if (reservoir) then - path%reservoir = trim(base_res) // "topology.data" - end if + if (reservoir) path%reservoir = trim(base_res) // "topology.data" ! Run prescan logic call prescan_inputs() @@ -77,12 +78,10 @@ subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) (nmax%inactive_residues == exp_inactive) if (.not. pass) then - print *, "Test FAILED" print *, " expected active =", exp_active, " got=", nmax%active_residues print *, " expected inactive =", exp_inactive, " got=", nmax%inactive_residues stop 1 - end if end subroutine run_test From 9b05eb7c5a0fe3b99694a763650dcce28c1554a6 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 11:23:08 +0100 Subject: [PATCH 11/13] minor improvement --- tests/readers/prescan/test_prescan_files.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/readers/prescan/test_prescan_files.f90 b/tests/readers/prescan/test_prescan_files.f90 index 54f805b..918b4fa 100644 --- a/tests/readers/prescan/test_prescan_files.f90 +++ b/tests/readers/prescan/test_prescan_files.f90 @@ -5,8 +5,8 @@ program test_prescan_files implicit none - character(len=LENPATH) :: base_main - character(len=LENPATH) :: base_res + character(len=LENPATH) :: base_main ! Path for main file + character(len=LENPATH) :: base_res ! Path for reservoir file !--------------------------------------------------------------------------- ! Test 1 : ZIF8-H2O without reservoir @@ -59,11 +59,11 @@ subroutine run_test(base_main, base_res, reservoir, exp_active, exp_inactive) integer, intent(in) :: exp_inactive ! Expected host residue ! Local variable - logical :: pass + logical :: pass ! Did the test pass ! Set paths - path%input = trim(base_main) // "input.maniac" - path%topology = trim(base_main) // "topology.data" + path%input = trim(base_main) // "input.maniac" + path%topology = trim(base_main) // "topology.data" path%parameters = trim(base_main) // "parameters.inc" status%reservoir_provided = reservoir From 638a1a39d78bd4f415a49fdf42401458db0216cd Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 11:40:24 +0100 Subject: [PATCH 12/13] Added tests for inputs --- .gitignore | 2 +- tests/readers/inputs/bad-input-01.maniac | 35 ++++++++--------- tests/readers/inputs/bad-input-02.maniac | 35 ++++++++--------- tests/readers/inputs/bad-input-03.maniac | 35 ++++++++--------- tests/readers/inputs/bad-input-04.maniac | 31 +++++++-------- tests/readers/inputs/good-input-01.maniac | 4 +- tests/readers/inputs/good-input-02.maniac | 46 +++++++++++++++++++++++ tests/readers/inputs/good-input-03.maniac | 44 ++++++++++++++++++++++ tests/readers/inputs/run-test.sh | 8 ++-- 9 files changed, 159 insertions(+), 81 deletions(-) create mode 100644 tests/readers/inputs/good-input-02.maniac create mode 100644 tests/readers/inputs/good-input-03.maniac diff --git a/.gitignore b/.gitignore index 971a84a..0fdbb0e 100644 --- a/.gitignore +++ b/.gitignore @@ -10,7 +10,7 @@ bin/* # Ignore produced files **/outputs/* **/*.dat -**/log.maniac +**/log*maniac # Ignore backup and temporary files (optional but useful) *~ diff --git a/tests/readers/inputs/bad-input-01.maniac b/tests/readers/inputs/bad-input-01.maniac index 5a7e50b..38b6236 100644 --- a/tests/readers/inputs/bad-input-01.maniac +++ b/tests/readers/inputs/bad-input-01.maniac @@ -1,23 +1,20 @@ -# Input for maniac program +# Input for maniac program with extra space before "translation_proba" -# Generic parameters -nb_block 1 # number of MC block -nb_step 1 # Number of MC step -temperature 300 # Temperature in Kelvin +nb_block 1 +nb_step 1 +temperature 300 -# Electrostatic interactions -ewald_tolerance 1e-5 # Tolerange long-range ewald -real_space_cutoff 17 # Cut-off electrostatic +ewald_tolerance 1e-5 +real_space_cutoff 17 -# Monte carlo move -translation_step 1 # Translational displacement (angstrom) -rotation_step_angle 0.685 # Rotational move (radian) -recalibrate_moves true # Enable automatic recalibration of move steps (true/false) +translation_step 1 +rotation_step_angle 0.685 +recalibrate_moves true - translation_proba 0.4 # Probability of attempting a translation move -rotation_proba 0.4 # Probability of attempting a rotation move -big_move_proba 0.0 # Probability of attempting a big move -insertion_deletion_proba 0.2 # Probability of attempting an insertion or a deletion + translation_proba 0.4 # Space before "translation_proba" +rotation_proba 0.4 +big_move_proba 0.0 +insertion_deletion_proba 0.2 # Define residues composition and activity begin_residue @@ -31,7 +28,7 @@ end_residue begin_residue name wat state actif - fugacity 50 + chemical_potential -4.5 types 12 13 14 names OW HW MW nb-atoms 4 @@ -40,8 +37,8 @@ end_residue begin_residue name moh state actif - fugacity 50 + chemical_potential -6.0 types 8 9 10 11 names HO HC CH OH nb-atoms 6 -end_residue \ No newline at end of file +end_residue diff --git a/tests/readers/inputs/bad-input-02.maniac b/tests/readers/inputs/bad-input-02.maniac index 58650aa..db7906d 100644 --- a/tests/readers/inputs/bad-input-02.maniac +++ b/tests/readers/inputs/bad-input-02.maniac @@ -1,23 +1,20 @@ -# Input for maniac program +# Input for maniac program with missing nb_block -# Generic parameters +# Missing nb_block +nb_step 1 +temperature 300 -nb_step 1 # Number of MC step -temperature 300 # Temperature in Kelvin +ewald_tolerance 1e-5 +real_space_cutoff 17 -# Electrostatic interactions -ewald_tolerance 1e-5 # Tolerange long-range ewald -real_space_cutoff 17 # Cut-off electrostatic +translation_step 1 +rotation_step_angle 0.685 +recalibrate_moves true -# Monte carlo move -translation_step 1 # Translational displacement (angstrom) -rotation_step_angle 0.685 # Rotational move (radian) -recalibrate_moves true # Enable automatic recalibration of move steps (true/false) - - translation_proba 0.4 # Probability of attempting a translation move -rotation_proba 0.4 # Probability of attempting a rotation move -big_move_proba 0.0 # Probability of attempting a big move -insertion_deletion_proba 0.2 # Probability of attempting an insertion or a deletion +translation_proba 0.4 +rotation_proba 0.4 +big_move_proba 0.0 +insertion_deletion_proba 0.2 # Define residues composition and activity begin_residue @@ -31,7 +28,7 @@ end_residue begin_residue name wat state actif - fugacity 50 + chemical_potential -4.5 types 12 13 14 names OW HW MW nb-atoms 4 @@ -40,8 +37,8 @@ end_residue begin_residue name moh state actif - fugacity 50 + chemical_potential -6.0 types 8 9 10 11 names HO HC CH OH nb-atoms 6 -end_residue \ No newline at end of file +end_residue diff --git a/tests/readers/inputs/bad-input-03.maniac b/tests/readers/inputs/bad-input-03.maniac index a544378..b47924a 100644 --- a/tests/readers/inputs/bad-input-03.maniac +++ b/tests/readers/inputs/bad-input-03.maniac @@ -1,23 +1,20 @@ -# Input for maniac program +# Input for maniac program with negative real_space_cutoff -# Generic parameters -nb_block 1 # number of MC block -nb_step 1 # Number of MC step -temperature 300 # Temperature in Kelvin +nb_block 1 +nb_step 1 +temperature 300 -# Electrostatic interactions -ewald_tolerance 1e-5 # Tolerange long-range ewald -real_space_cutoff -17 # Cut-off electrostatic +ewald_tolerance 1e-5 +real_space_cutoff -17 -# Monte carlo move -translation_step 1 # Translational displacement (angstrom) -rotation_step_angle 0.685 # Rotational move (radian) -recalibrate_moves true # Enable automatic recalibration of move steps (true/false) +translation_step 1 +rotation_step_angle 0.685 +recalibrate_moves true - translation_proba 0.4 # Probability of attempting a translation move -rotation_proba 0.4 # Probability of attempting a rotation move -big_move_proba 0.0 # Probability of attempting a big move -insertion_deletion_proba 0.2 # Probability of attempting an insertion or a deletion +translation_proba 0.4 +rotation_proba 0.4 +big_move_proba 0.0 +insertion_deletion_proba 0.2 # Define residues composition and activity begin_residue @@ -31,7 +28,7 @@ end_residue begin_residue name wat state actif - fugacity 50 + chemical_potential -4.5 types 12 13 14 names OW HW MW nb-atoms 4 @@ -40,8 +37,8 @@ end_residue begin_residue name moh state actif - fugacity 50 + chemical_potential -6.0 types 8 9 10 11 names HO HC CH OH nb-atoms 6 -end_residue \ No newline at end of file +end_residue diff --git a/tests/readers/inputs/bad-input-04.maniac b/tests/readers/inputs/bad-input-04.maniac index fe888ee..f2a1e2a 100644 --- a/tests/readers/inputs/bad-input-04.maniac +++ b/tests/readers/inputs/bad-input-04.maniac @@ -1,23 +1,20 @@ -# Input for maniac program +# Input for maniac program with missing chemical potential -# Generic parameters -nb_block 1 # number of MC block -nb_step 1 # Number of MC step -temperature 300 # Temperature in Kelvin +nb_block 1 +nb_step 1 +temperature 300 -# Electrostatic interactions -ewald_tolerance 1e-5 # Tolerange long-range ewald -real_space_cutoff 17 # Cut-off electrostatic +ewald_tolerance 1e-5 +real_space_cutoff 17 -# Monte carlo move -translation_step 1 # Translational displacement (angstrom) -rotation_step_angle 0.685 # Rotational move (radian) -recalibrate_moves true # Enable automatic recalibration of move steps (true/false) +translation_step 1 +rotation_step_angle 0.685 +recalibrate_moves true -translation_proba 0.4 # Probability of attempting a translation move -rotation_proba 0.4 # Probability of attempting a rotation move -big_move_proba 0.0 # Probability of attempting a big move -insertion_deletion_proba 0.2 # Probability of attempting an insertion or a deletion +translation_proba 0.4 +rotation_proba 0.4 +big_move_proba 0.0 +insertion_deletion_proba 0.2 # Define residues composition and activity begin_residue @@ -31,7 +28,7 @@ end_residue begin_residue name wat state actif - fugacity 50 + chemical_potential -4.5 types 12 13 14 names OW HW MW nb-atoms 4 diff --git a/tests/readers/inputs/good-input-01.maniac b/tests/readers/inputs/good-input-01.maniac index 93f6dee..0739b19 100644 --- a/tests/readers/inputs/good-input-01.maniac +++ b/tests/readers/inputs/good-input-01.maniac @@ -30,7 +30,7 @@ end_residue begin_residue name wat state actif - fugacity 50 + chemical_potential -4.5 types 12 13 14 names OW HW MW nb-atoms 4 @@ -39,7 +39,7 @@ end_residue begin_residue name moh state actif - fugacity 50 + chemical_potential -6.0 types 8 9 10 11 names HO HC CH OH nb-atoms 6 diff --git a/tests/readers/inputs/good-input-02.maniac b/tests/readers/inputs/good-input-02.maniac new file mode 100644 index 0000000..e5b4af1 --- /dev/null +++ b/tests/readers/inputs/good-input-02.maniac @@ -0,0 +1,46 @@ +# Input for maniac program + +# Generic parameters +nb_block 1 +nb_step 1 +temperature 300 + +# Electrostatic interactions +ewald_tolerance 1e-5 +real_space_cutoff 17 + +# Monte carlo move +translation_step 1 +rotation_step_angle 0.685 +recalibrate_moves true + +translation_proba 0.4 +rotation_proba 0.4 +widom_proba 0.2 + +# Define residues composition and activity +begin_residue + name zif + state inactif + types 1 2 3 4 5 6 7 + names C1 C2 C3 H2 H3 N Zn + nb-atoms 2208 +end_residue + +begin_residue + name wat + state actif + chemical_potential -4.5 + types 12 13 14 + names OW HW MW + nb-atoms 4 +end_residue + +begin_residue + name moh + state actif + chemical_potential -6.0 + types 8 9 10 11 + names HO HC CH OH + nb-atoms 6 +end_residue diff --git a/tests/readers/inputs/good-input-03.maniac b/tests/readers/inputs/good-input-03.maniac new file mode 100644 index 0000000..361fb30 --- /dev/null +++ b/tests/readers/inputs/good-input-03.maniac @@ -0,0 +1,44 @@ +# Input for maniac program + +# Generic parameters +nb_block 1 +nb_step 1 +temperature 300 + +# Electrostatic interactions +ewald_tolerance 1e-5 +real_space_cutoff 17 + +# Monte carlo move +translation_step 1 +rotation_step_angle 0.685 +recalibrate_moves true + +translation_proba 0.4 +rotation_proba 0.4 +widom_proba 0.2 + +# Define residues composition and activity +begin_residue + name zif + state inactif + types 1 2 3 4 5 6 7 + names C1 C2 C3 H2 H3 N Zn + nb-atoms 2208 +end_residue + +begin_residue + name wat + state actif + types 12 13 14 + names OW HW MW + nb-atoms 4 +end_residue + +begin_residue + name moh + state actif + types 8 9 10 11 + names HO HC CH OH + nb-atoms 6 +end_residue diff --git a/tests/readers/inputs/run-test.sh b/tests/readers/inputs/run-test.sh index 36df301..e03c5d2 100755 --- a/tests/readers/inputs/run-test.sh +++ b/tests/readers/inputs/run-test.sh @@ -15,10 +15,10 @@ for input in good-input-*.maniac; do mkdir -p "$outputs" # Run silently (everything goes into log.maniac) - $build_path -i "$input" -d "$data" -p "$inc" -o "$outputs" > log.maniac 2>&1 + $build_path -i "$input" -d "$data" -p "$inc" -o "$outputs" > "log.$input" 2>&1 # === Test: program termination - if grep -q "MANIAC-MC simulation completed " log.maniac; then + if grep -q "MANIAC-MC simulation completed " "log.$input"; then echo "✅ [PASS] $input : Simulation terminated normally" else echo "❌ [FAIL] $input : No termination message found" @@ -32,9 +32,9 @@ for input in bad-input-*.maniac; do rm -rf "$outputs" mkdir -p "$outputs" - $build_path -i "$input" -d "$data" -p "$inc" -o "$outputs" > log.maniac 2>&1 || true # don't exit on error + $build_path -i "$input" -d "$data" -p "$inc" -o "$outputs" > "log.$input" 2>&1 || true # don't exit on error - if grep -E -q "(Error|STOP)" log.maniac; then + if grep -E -q "(Error|STOP)" "log.$input"; then echo "✅ [PASS] $input : Correctly failed with error" else echo "❌ [FAIL] $input : Unexpectedly succeeded (no error in log)" From b59df747f8fc2f0810ca61f3250aa449ead60e15 Mon Sep 17 00:00:00 2001 From: Simon Gravelle Date: Sun, 30 Nov 2025 12:11:18 +0100 Subject: [PATCH 13/13] fixed troncated error --- src/prepare_utils.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/prepare_utils.f90 b/src/prepare_utils.f90 index 43d8e40..77278d6 100644 --- a/src/prepare_utils.f90 +++ b/src/prepare_utils.f90 @@ -61,8 +61,10 @@ subroutine allocate_array() ! Allocate complex arrays for wave vector components kmax_max = maxval(ewald%param%kmax) - allocate(ewald%phase%factor_host(3, res%number, 0:nmax%inactive_residues, 1:nmax%atoms_inactive_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? - allocate(ewald%phase%factor_guest(3, res%number, 0:nmax%active_residues, 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) ! TOFIX, do not use NB_MAX_MOLECULE systematical ? + allocate(ewald%phase%factor_host(3, res%number, 0:nmax%inactive_residues, & + 1:nmax%atoms_inactive_residue, -kmax_max:kmax_max)) + allocate(ewald%phase%factor_guest(3, res%number, 0:nmax%active_residues, & + 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) allocate(ewald%phase%factor_old(3, 1:nmax%atoms_active_residue, -kmax_max:kmax_max)) ! Allocate temporary arrays once