Consistent use of chemical potential (#24)

* Started changing units for chemical potential

* progressed in using chemical potential

* keep improving units management

* Improved compute_acceptance_probability

* updated MC

* Updated Widom

* improved move writer

* Improved writting

* Improved and cleaned constant

* prepare releavse of 0.3.3-beta

Author: simongravelle

mc-topology

@@ -1,16 +1,25 @@
 module constants
 
     use, intrinsic :: iso_fortran_env, only: real64
+    use parameters
 
     implicit none
 
-    real(real64), parameter :: PI = 3.14159265358979323846_real64 ! Constant Pi
-    real(real64), parameter :: TWOPI = 2.0_real64 * PI ! Two times Pi
-    real(real64), parameter :: SQRTPI = SQRT(PI) ! Square root of Pi
-    real(real64), parameter :: KB_JK = 1.380658D-23 ! Boltzmann constant (J/K)
-    real(real64), parameter :: KB_kcalmol = 0.0019872041_real64 ! Boltzmann constant (kcal/mol)
-    real(real64), parameter :: EPS0_INV_kcalA = 332.0637_real64 ! e^2 / 4 pi epsilon_0 (kcal Å/(mol e^2))
-    real(real64), parameter :: KB_eVK = 8.6173852D-5    ! Boltzmann constant (eV/K)
+    real(real64), parameter :: PI = 3.1415926536_real64     ! Constant Pi
+    real(real64), parameter :: TWOPI = 2.0_real64 * PI      ! Two times Pi
+    real(real64), parameter :: SQRTPI = SQRT(PI)            ! Square root of Pi
+    real(real64), parameter :: H_PLANCK = 6.62607015d-34    ! Planck constant (J s)
+    real(real64), parameter :: EPS0 = 8.854187817d-12       ! F/m = C^2/(N m^2)
+    real(real64), parameter :: NA = 6.02214076e23           ! Na (mol-1)
+    real(real64), parameter :: KB = 1.380658d-23            ! Boltzmann constant (J/K)
+    real(real64), parameter :: E_CHARGE  = 1.602176634d-19  ! Elementary charge (C)
+
+    real(real64), parameter :: EPS0_INV = E_CHARGE**2 / (4.0_real64 * PI * EPS0) ! Coulomb prefactor e^2 / 4 pi epsilon_0, SI units, N m^2
+    real(real64), parameter :: EPS0_INV_real = EPS0_INV * J_to_kcal * m_to_A * NA ! Coulomb prefactor in real units, 
+
+    real(real64) :: beta                                    ! 1 / (kB T)
+    real(real64), parameter :: KB_kcalmol = KB * NA * J_to_kcal ! Boltzmann constant (in kcal/mol)
+
     real(real64), parameter :: zero = 0.0_real64        ! Zero in real64
     real(real64), parameter :: quarter = 0.25_real64    ! 1/4 in real64
     real(real64), parameter :: half = 0.5_real64        ! 1/2 in real64
@@ -19,8 +28,5 @@ module constants
     real(real64), parameter :: three = 3.0_real64       ! Three in real64
     real(real64), parameter :: four = 4.0_real64        ! Four in real64
     real(real64), parameter :: ten = 10.0_real64        ! Ten in real64
-    real(real64), parameter :: error = 1.0D-10          ! Small number for error calculation
-    real(real64), parameter :: Hplank = 6.62607015D-34  ! Planck constant (J s)
-    real(real64), parameter :: HBAR = 1.05457182D-34    ! hbar (m^2 kg / s)
-    real(real64), parameter :: NA = 6.022e23            ! Na (mol-1)
+
 end module constants

src/constants.f90

@@ -1461,7 +1461,7 @@ subroutine TransformCoordinate(box)
                                     com(1), com(2), com(3))
 
                     total_mass = ComputeMass(nb%atom_in_residue(i), tmp_atom_masses_1d)
-                    res%mass_residue(i) = total_mass 
+                    res%mass(i) = total_mass 
 
                     ! Save original CoM before applying periodic boundary conditions
                     original_com = com ! Save the CoM

src/data_parser.f90

@@ -189,7 +189,7 @@ subroutine SingleMolPairwiseEnergy(box, residue_type_1, molecule_index_1, e_non_
         end do
 
         ! Re-scale energy
-        e_coulomb = e_coulomb * EPS0_INV_kcalA   ! In kcal/mol
+        e_coulomb = e_coulomb * EPS0_INV_real   ! In kcal/mol
 
     end subroutine SingleMolPairwiseEnergy
 
@@ -379,7 +379,7 @@ subroutine SingleMolEwaldSelf(residue_type, self_energy_1)
         end do
 
         ! Convert to kcal/mol at the end
-        self_energy_1 = self_energy_1 * EPS0_INV_kcalA  ! In kcal/mol
+        self_energy_1 = self_energy_1 * EPS0_INV_real  ! In kcal/mol
 
     end subroutine SingleMolEwaldSelf
 
@@ -455,7 +455,7 @@ subroutine ComputePairInteractionEnergy_singlemol(box, residue_type_1, molecule_
         end do
 
         ! Convert to kcal/mol at the end
-        e_coulomb = e_coulomb * EPS0_INV_kcalA                          ! In kcal/mol
+        e_coulomb = e_coulomb * EPS0_INV_real                          ! In kcal/mol
 
     end subroutine ComputePairInteractionEnergy_singlemol
 

src/energy_utils.f90

@@ -132,8 +132,7 @@ subroutine ComputeReciprocalEnergy(u_recipCoulomb)
         end do
 
         ! Convert accumulated energy to correct units (kcal/mol)
-        ! e^2 x Å^2 * kcal Å / (mol e^2 Å^3) = kcal/mol
-        u_recipCoulomb = u_recipCoulomb * EPS0_INV_kcalA * TWOPI / primary%volume ! In kcal/mol
+        u_recipCoulomb = u_recipCoulomb * EPS0_INV_real * TWOPI / primary%volume ! In kcal/mol
 
     end subroutine ComputeReciprocalEnergy
 
@@ -253,10 +252,8 @@ subroutine ComputeRecipEnergySingleMol(residue_type, molecule_index, u_recipCoul
 
         end do
 
-        !----------------------------------------------
         ! Convert accumulated energy to physical units:
-        !----------------------------------------------
-        u_recipCoulomb_new = u_recipCoulomb_new * EPS0_INV_kcalA * TWOPI / primary%volume ! In kcal/mol
+        u_recipCoulomb_new = u_recipCoulomb_new * EPS0_INV_real * TWOPI / primary%volume ! In kcal/mol
 
     end subroutine ComputeRecipEnergySingleMol
 
@@ -306,7 +303,7 @@ subroutine ComputeEwaldSelfInteractionSingleMol(residue_type, self_energy)
             self_energy = self_energy - ewald%alpha / SQRTPI * charge_1**2
         end do
 
-        self_energy = self_energy * EPS0_INV_kcalA ! In kcal/mol
+        self_energy = self_energy * EPS0_INV_real ! In kcal/mol
 
     end subroutine ComputeEwaldSelfInteractionSingleMol
 
@@ -366,7 +363,7 @@ subroutine ComputeIntraResidueRealCoulombEnergySingleMol(residue_type, molecule_
             end do
         end do
 
-        u_intraCoulomb = u_intraCoulomb * EPS0_INV_kcalA ! In kcal/mol
+        u_intraCoulomb = u_intraCoulomb * EPS0_INV_real ! In kcal/mol
 
     end subroutine ComputeIntraResidueRealCoulombEnergySingleMol
 

src/ewald_energy.f90

@@ -175,4 +175,42 @@ pure real(real64) function vector_norm(v)
 
     end function vector_norm
 
+    !--------------------------------------------------------------------
+    ! Pads the given column title to 12 characters and adds one space (A12,1X),
+    ! ensuring consistent alignment with numeric columns printed using I12.
+    ! Appends the formatted field to the growing header line.
+    !--------------------------------------------------------------------
+    subroutine add_col(line, title)
+
+        ! Input argument
+        character(len=*), intent(inout) :: line
+        character(len=*), intent(in)    :: title
+
+        ! Local variable
+        character(len=16) :: tmp
+
+        ! EXACT match to "(1X,I12)" used by data columns
+        write(tmp,'(1X, A12)') title
+
+        ! DO NOT trim, append raw padded field
+        line = line // tmp
+
+    end subroutine add_col
+
+    subroutine add_first_col(line, title)
+
+        ! Input arguments
+        character(len=*), intent(inout) :: line
+        character(len=*), intent(in)    :: title
+
+        ! Local variable
+        character(len=16) :: tmp
+
+        ! FIRST column: matches (I12)
+        write(tmp,'(A12)') title
+        line = line // tmp
+
+    end subroutine add_first_col
+
+
 end module helper_utils

src/helper_utils.f90

@@ -259,7 +259,8 @@ subroutine AllocateAtomArrays()
         allocate(reservoir%mol_com(3, nb%type_residue, NB_MAX_MOLECULE))
         allocate(reservoir%site_offset(3, nb%type_residue, NB_MAX_MOLECULE, nb%max_atom_in_residue))
         allocate(reservoir%num_residues(nb%type_residue))
-        allocate(res%mass_residue(nb%type_residue))
+        allocate(res%mass(nb%type_residue))
+        allocate(res%lambda(nb%type_residue))
 
         ! Allocate parameter arrays
         allocate(coeff%sigma(nb%type_residue, nb%type_residue,&
@@ -270,6 +271,7 @@ subroutine AllocateAtomArrays()
         allocate(res%types_2d(nb%type_residue, nb%max_type_per_residue))
         allocate(res%names_2d(nb%type_residue, nb%max_type_per_residue))
         allocate(input%fugacity(nb%type_residue))
+        allocate(input%chemical_potential(nb%type_residue))
         allocate(nb%types_per_residue(nb%type_residue))
         allocate(nb%atom_in_residue(nb%type_residue))
         allocate(nb%bonds_per_residue(nb%type_residue))
@@ -325,6 +327,7 @@ subroutine ParseInputFile(INFILE)
         logical :: has_translation_step, has_rotation_step, has_recalibrate
         logical :: has_translation_proba, has_rotation_proba
         logical :: has_insertdel_proba, has_swap_proba, has_widom_proba
+        logical :: has_fugacity, has_chemical_potential
 
         has_nb_block = .false.
         has_nb_step = .false.
@@ -347,6 +350,8 @@ subroutine ParseInputFile(INFILE)
 
         ! Initialize all fugacities to -1.0 (indicating unset)
         input%fugacity(:) = -one
+        ! Initialize all chemical_potential to 0.0 (indicating unset)
+        input%chemical_potential(:) = zero
 
         do
             read(INFILE, '(A)', IOSTAT=ios) line
@@ -379,7 +384,7 @@ subroutine ParseInputFile(INFILE)
                 read(rest_line, *, iostat=ios) val_real
                 if (ios /= 0) error stop "Error reading temperature"
                 if (val_real <= zero) error stop "Invalid temperature: must be > 0"
-                input%temp_K = val_real
+                input%temperature = val_real
                 has_temp = .true.
 
             case ("seed")
@@ -473,7 +478,7 @@ subroutine ParseInputFile(INFILE)
             end select
 
             ! ===============================
-            ! Residue parsing (types, names, fugacity, nb-atoms)
+            ! Residue parsing (types, names, fugacity, chemical potential nb-atoms)
             ! ===============================
             if (in_residue_block) then
 
@@ -520,6 +525,11 @@ subroutine ParseInputFile(INFILE)
                     read(rest_line, *, iostat=ios) val_real
                     input%fugacity(nb%type_residue + 1) = val_real
 
+                ! Check if the line specifies the chemical potential
+                else if (trim(token) == 'chemical_potential') then
+                    read(rest_line, *, iostat=ios) val_real
+                    input%chemical_potential(nb%type_residue + 1) = val_real
+
                 ! Check if the line specifies the number of atoms
                 else if (trim(token) == 'nb-atoms') then
                     read(rest_line, *, iostat=ios) val_int
@@ -571,12 +581,24 @@ subroutine ParseInputFile(INFILE)
 
         ! Validate fugacity for active residues
         do val_int = 1, nb%type_residue
-            if (input%is_active(val_int) == 1) then
-                if (input%fugacity(val_int) < zero) then
-                    call AbortRun("Fugacity not provided or invalid for active residue: " &
-                        // trim(res%names_1d(val_int)))
-                end if
+
+            if (input%is_active(val_int) == 0) cycle
+
+            has_fugacity = (input%fugacity(val_int) >= zero)
+            has_chemical_potential = (input%chemical_potential(val_int) < zero)
+
+            ! Rule 1: at least one must be provided
+            if (.not.(has_fugacity .or. has_chemical_potential)) then
+                call AbortRun("Neither fugacity nor chemical potential provided for active residue: " // &
+                            trim(res%names_1d(val_int)))
+            end if
+
+            ! Rule 2: cannot both be provided
+            if (has_fugacity .and. has_chemical_potential) then
+                call AbortRun("Both fugacity and chemical potential were specified for active residue: " // &
+                            trim(res%names_1d(val_int)))
             end if
+
         end do
 
         ! === Validation of required parameters ===
@@ -628,7 +650,7 @@ subroutine SortResidues()
         integer, allocatable :: keys(:), order(:)
         character(len=10), allocatable :: tmp_names_1d(:)
         integer, allocatable :: tmp_is_active(:), tmp_atom_in_residue(:), tmp_types_per_residue(:)
-        real(real64), allocatable :: tmp_fugacity(:)
+        real(real64), allocatable :: tmp_fugacity(:), tmp_chemical_potential(:)
         integer, allocatable :: tmp_types_2d(:,:)
         character(len=10), allocatable :: tmp_names_2d(:,:)
 
@@ -640,6 +662,7 @@ subroutine SortResidues()
         allocate(tmp_names_1d(n))
         allocate(tmp_is_active(n))
         allocate(tmp_fugacity(n))
+        allocate(tmp_chemical_potential(n))
         allocate(tmp_atom_in_residue(n))
         allocate(tmp_types_per_residue(n))
         allocate(tmp_types_2d(size(res%types_2d,1), size(res%types_2d,2)))
@@ -670,6 +693,7 @@ subroutine SortResidues()
             tmp_names_1d(i)          = res%names_1d(k)
             tmp_is_active(i)         = input%is_active(k)
             tmp_fugacity(i)          = input%fugacity(k)
+            tmp_chemical_potential(i) = input%chemical_potential(k)
             tmp_atom_in_residue(i)   = nb%atom_in_residue(k)
             tmp_types_per_residue(i) = nb%types_per_residue(k)
             tmp_types_2d(i,:)        = res%types_2d(k,:)
@@ -680,13 +704,14 @@ subroutine SortResidues()
         res%names_1d          = tmp_names_1d
         input%is_active       = tmp_is_active
         input%fugacity        = tmp_fugacity
+        input%chemical_potential = tmp_chemical_potential
         nb%atom_in_residue    = tmp_atom_in_residue
         nb%types_per_residue  = tmp_types_per_residue
         res%types_2d          = tmp_types_2d
         res%names_2d          = tmp_names_2d
 
         ! Deallocate temporary arrays
-        deallocate(keys, order, tmp_names_1d, tmp_is_active, tmp_fugacity, &
+        deallocate(keys, order, tmp_names_1d, tmp_is_active, tmp_fugacity, tmp_chemical_potential, &
                 tmp_atom_in_residue, tmp_types_per_residue, tmp_types_2d, tmp_names_2d)
 
     end subroutine SortResidues