Commit 17bc731c authored by Remco Hens's avatar Remco Hens
Browse files

fix bug enthalpy insertion

parent 2deef8c2
No preview for this file type
......@@ -14,14 +14,16 @@ ionic systems). For documentation, we would like to refer to the following:
- [Main publication of Brick-CFCMC][paper1] and its Supporting Information (open
access)
- [Review paper on the CFCMC technique][paper2] (open access)
- [Thermodynamic Integrations and Collective Trial Moves in Brick-CFCMC][paper3] (open access)
- [PhD thesis of Remco Hens][thesis1] (open access) [(alternative website)][thesis1a]
- [PhD thesis of Ahmadreza Rahbari][thesis2] (open access) [(alternative website)][thesis2a]
-[Website][web1] of the Engineering Thermodynamics group at Delft University of Technology
- [Website][web1] of the Engineering Thermodynamics group at Delft University of Technology
[gitlabeth]: <https://gitlab.com/ETh_TU_Delft/Brick-CFCMC>
[ethmanual]: <https://homepage.tudelft.nl/v9k6y/Brick-CFCMC/Brick-CFCMC.pdf>
[paper1]: <https://doi.org/10.1021/acs.jcim.0c00334>
[paper2]: <https://doi.org/10.1080/08927022.2020.1828585>
[paper3]: <https://doi.org/10.1021/acs.jcim.1c00652>
[thesis1]: <https://doi.org/10.4233/uuid:41c32a8f-2db7-4091-abb5-4d6a5e596345>
[thesis1a]: <https://homepage.tudelft.nl/v9k6y/thesis-Hens.pdf>
[thesis2]: <https://doi.org/10.4233/uuid:eb04d860-281a-4c6b-8c5b-263f526d0bd9>
......@@ -32,6 +34,7 @@ If you use Brick-CFCMC in your publication, please cite the following papers:
- [J. Chem. Inf. Model., 2020, 60, 2678-2682][paper1]
- [Molecular Simulation, 2021, 47, 804-823][paper2]
- [J. Chem. Inf. Model., 2021, 61, 3752-3757][paper3]
# Quick Start
We assume that you are working on a Linux/Unix system, and that you
......
SUBROUTINE Enthalpy_insertion(Ichoice)
implicit none
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC Test particle insertion method by Ciccotti and Frenkel CCC
CCC for partial molar enthalpy and partial molar volume CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
include "global_variables.inc"
include "output.inc"
integer Ninsertions
Parameter (Ninsertions = 10)
integer Ichoice,I,J,Ib,Tm,Imol
double precision E_LJ_Inter,E_LJ_Intra,E_EL_Real,E_EL_Intra,E_Bending,E_Torsion,E_EL_Excl,dU,U
double precision dX,dY,dZ,Ran_Uniform
double precision H_test_1(2,MaxMolType),H_test_2(2,MaxMolType),H_test_3(2,MaxMolType),
& V_test_1(2,MaxMolType),V_test_3(2,MaxMolType),Counter,
& H2_test_1(2,MaxMolType),H2_test_3(2,MaxMolType)
double precision dU_dl
logical L_Overlap_Inter,L_Overlap_Intra
save H_test_1,H_test_2,H_test_3,V_test_1,V_test_3,H2_test_1,H2_test_3,Counter
CCC INITIALIZE
IF(Ichoice.EQ.0) THEN
DO Ib=1,N_Box
DO Tm=1,N_MolType
H_test_1(Ib,Tm) = 0.0d0
H_test_2(Ib,Tm) = 0.0d0
H_test_3(Ib,Tm) = 0.0d0
V_test_1(Ib,Tm) = 0.0d0
V_test_3(Ib,Tm) = 0.0d0
H2_test_1(Ib,Tm) = 0.0d0
H2_test_3(Ib,Tm) = 0.0d0
END DO
END DO
Counter = 0.0d0
CCC SAMPLE
ELSEIF(Ichoice.EQ.1) THEN
DO J=1,Ninsertions
DO Ib=1,N_Box
DO Tm=1,N_MolType
Imol = N_MolInBox(Ib) + 1
Ibox(Imol) = Ib
TypeMol(Imol) = Tm
L_Frac(Imol) = .false.
XCM(Imol) = 0.0d0
YCM(Imol) = 0.0d0
ZCM(Imol) = 0.0d0
dX = Ran_Uniform()*BoxSize(Ib)
dY = Ran_Uniform()*BoxSize(Ib)
dZ = Ran_Uniform()*BoxSize(Ib)
DO I=1,N_AtomInMolType(Tm)
X(Imol,I) = Xin(Tm,I) + dX
Y(Imol,I) = Yin(Tm,I) + dY
Z(Imol,I) = Zin(Tm,I) + dZ
XCM(Imol) = XCM(Imol) + X(Imol,I)
YCM(Imol) = YCM(Imol) + Y(Imol,I)
ZCM(Imol) = ZCM(Imol) + Z(Imol,I)
END DO
XCM(Imol) = XCM(Imol)/dble(N_AtomInMolType(Tm))
YCM(Imol) = YCM(Imol)/dble(N_AtomInMolType(Tm))
ZCM(Imol) = ZCM(Imol)/dble(N_AtomInMolType(Tm))
CALL Random_Orientation(Imol)
CALL Place_molecule_back_in_box(Imol)
CALL Energy_Molecule(Imol,E_LJ_Inter,E_LJ_Intra,E_EL_Real,E_EL_Intra,E_EL_Excl,
& E_Bending,E_Torsion,L_Overlap_Inter,L_Overlap_Intra,dU_dl)
IF(.NOT.(L_Overlap_Inter.OR.L_Overlap_Intra)) THEN
U = U_LJ_Inter(Ib) + U_LJ_Tail(Ib)
& + U_EL_Real(Ib) + U_EL_Excl(Ib) + U_EL_Self(Ib) + U_EL_Four(Ib)
& + U_LJ_Intra(Ib) + U_EL_Intra(Ib) + U_Bending_Total(Ib) + U_Torsion_Total(Ib)
dU=E_LJ_Inter+E_LJ_Intra+E_EL_Real+E_EL_Intra+E_EL_Excl+E_Bending+E_Torsion
H_test_1(Ib,Tm) = H_test_1(Ib,Tm)
& + (U + dU + Pressure*Volume(Ib))*Volume(Ib)*dexp(-beta*dU)
H_test_2(Ib,Tm) = H_test_2(Ib,Tm) + Volume(Ib)*dexp(-beta*dU)
V_test_1(Ib,Tm) = V_test_1(Ib,Tm) + Volume(Ib)*Volume(Ib)*dexp(-beta*dU)
H2_test_1(Ib,Tm) = H2_test_1(Ib,Tm) + (U + dU + Pressure*Volume(Ib))*
& (U + dU + Pressure*Volume(Ib))*Volume(Ib)*dexp(-beta*dU)
END IF
H_test_3(Ib,Tm) = H_test_3(Ib,Tm) + (U + Pressure*Volume(Ib))
V_test_3(Ib,Tm) = V_test_3(Ib,Tm) + Volume(Ib)
H2_test_3(Ib,Tm) = H2_test_3(Ib,Tm)
& + (U + Pressure*Volume(Ib))*(U + Pressure*Volume(Ib))
END DO
END DO
END DO
Counter = Counter + dble(Ninsertions)
CCC WRITE
ELSEIF(Ichoice.EQ.2) THEN
WRITE(6,'(A)') "Enthalpy"
DO Tm=1,N_MolType
WRITE(6,'(1x,A19,18x,2e20.10e3)') C_MolType(Tm),
& (-1.0d0/beta + H_test_1(Ib,Tm)/H_test_2(Ib,Tm) - H_test_3(Ib,Tm)/Counter, Ib=1,N_Box)
END DO
WRITE(6,*)
WRITE(6,'(A)') "Partial volume"
DO Tm=1,N_MolType
WRITE(6,'(1x,A19,18x,2e20.10e3)') C_MolType(Tm),
& (V_test_1(Ib,Tm)/H_test_2(Ib,Tm) - V_test_3(Ib,Tm)/Counter, Ib=1,N_Box)
END DO
WRITE(6,*)
WRITE(6,'(A)') "Heat Capacity"
DO Tm=1,N_MolType
WRITE(6,'(1x,A19,18x,2e20.10e3)') C_MolType(Tm),
& (1.0d0/(beta*beta) + (H_test_1(Ib,Tm)/H_test_2(Ib,Tm))*(H_test_1(Ib,Tm)/H_test_2(Ib,Tm))
& - H2_test_1(Ib,Tm)/H_test_2(Ib,Tm) - (H_test_3(Ib,Tm)/Counter)*(H_test_3(Ib,Tm)/Counter)
& + H2_test_3(Ib,Tm)/Counter , Ib=1,N_Box)
END DO
END IF
RETURN
END
SUBROUTINE Enthalpy_insertion(Ichoice)
implicit none
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC Test particle insertion method by Ciccotti and Frenkel CCC
CCC for partial molar enthalpy and partial molar volume CCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
include "global_variables.inc"
include "output.inc"
integer Ninsertions
Parameter (Ninsertions = 10)
integer Ichoice,I,J,Ib,Tm,Imol
double precision E_LJ_Inter,E_LJ_Intra,E_EL_Real,E_EL_Intra,E_Bending,E_Torsion,E_EL_Excl,dU,U
double precision dX,dY,dZ,Ran_Uniform
double precision H_test_1(2,MaxMolType),H_test_2(2,MaxMolType),H_test_3(2,MaxMolType),
& V_test_1(2,MaxMolType),V_test_3(2,MaxMolType),Counter,
& H2_test_1(2,MaxMolType),H2_test_3(2,MaxMolType)
double precision dU_dl
logical L_Overlap_Inter,L_Overlap_Intra
save H_test_1,H_test_2,H_test_3,V_test_1,V_test_3,H2_test_1,H2_test_3,Counter
CCC INITIALIZE
IF(Ichoice.EQ.0) THEN
DO Ib=1,N_Box
DO Tm=1,N_MolType
H_test_1(Ib,Tm) = 0.0d0
H_test_2(Ib,Tm) = 0.0d0
H_test_3(Ib,Tm) = 0.0d0
V_test_1(Ib,Tm) = 0.0d0
V_test_3(Ib,Tm) = 0.0d0
H2_test_1(Ib,Tm) = 0.0d0
H2_test_3(Ib,Tm) = 0.0d0
END DO
END DO
Counter = 0.0d0
CCC SAMPLE
ELSEIF(Ichoice.EQ.1) THEN
DO J=1,Ninsertions
DO Ib=1,N_Box
DO Tm=1,N_MolType
Imol = N_MolInBox(Ib) + 1
Ibox(Imol) = Ib
TypeMol(Imol) = Tm
L_Frac(Imol) = .false.
XCM(Imol) = 0.0d0
YCM(Imol) = 0.0d0
ZCM(Imol) = 0.0d0
dX = Ran_Uniform()*BoxSize(Ib)
dY = Ran_Uniform()*BoxSize(Ib)
dZ = Ran_Uniform()*BoxSize(Ib)
DO I=1,N_AtomInMolType(Tm)
X(Imol,I) = Xin(Tm,I) + dX
Y(Imol,I) = Yin(Tm,I) + dY
Z(Imol,I) = Zin(Tm,I) + dZ
XCM(Imol) = XCM(Imol) + X(Imol,I)
YCM(Imol) = YCM(Imol) + Y(Imol,I)
ZCM(Imol) = ZCM(Imol) + Z(Imol,I)
END DO
XCM(Imol) = XCM(Imol)/dble(N_AtomInMolType(Tm))
YCM(Imol) = YCM(Imol)/dble(N_AtomInMolType(Tm))
ZCM(Imol) = ZCM(Imol)/dble(N_AtomInMolType(Tm))
CALL Random_Orientation(Imol)
CALL Place_molecule_back_in_box(Imol)
CALL Energy_Molecule(Imol,E_LJ_Inter,E_LJ_Intra,E_EL_Real,E_EL_Intra,E_EL_Excl,
& E_Bending,E_Torsion,L_Overlap_Inter,L_Overlap_Intra,dU_dl)
U = U_LJ_Inter(Ib) + U_LJ_Tail(Ib)
& + U_EL_Real(Ib) + U_EL_Excl(Ib) + U_EL_Self(Ib) + U_EL_Four(Ib)
& + U_LJ_Intra(Ib) + U_EL_Intra(Ib) + U_Bending_Total(Ib) + U_Torsion_Total(Ib)
IF(.NOT.(L_Overlap_Inter.OR.L_Overlap_Intra)) THEN
dU=E_LJ_Inter+E_LJ_Intra+E_EL_Real+E_EL_Intra+E_EL_Excl+E_Bending+E_Torsion
H_test_1(Ib,Tm) = H_test_1(Ib,Tm)
& + (U + dU + Pressure*Volume(Ib))*Volume(Ib)*dexp(-beta*dU)
H_test_2(Ib,Tm) = H_test_2(Ib,Tm) + Volume(Ib)*dexp(-beta*dU)
V_test_1(Ib,Tm) = V_test_1(Ib,Tm) + Volume(Ib)*Volume(Ib)*dexp(-beta*dU)
H2_test_1(Ib,Tm) = H2_test_1(Ib,Tm) + (U + dU + Pressure*Volume(Ib))*
& (U + dU + Pressure*Volume(Ib))*Volume(Ib)*dexp(-beta*dU)
END IF
H_test_3(Ib,Tm) = H_test_3(Ib,Tm) + (U + Pressure*Volume(Ib))
V_test_3(Ib,Tm) = V_test_3(Ib,Tm) + Volume(Ib)
H2_test_3(Ib,Tm) = H2_test_3(Ib,Tm)
& + (U + Pressure*Volume(Ib))*(U + Pressure*Volume(Ib))
END DO
END DO
END DO
Counter = Counter + dble(Ninsertions)
CCC WRITE
ELSEIF(Ichoice.EQ.2) THEN
WRITE(6,'(A)') "Enthalpy"
DO Tm=1,N_MolType
WRITE(6,'(1x,A19,18x,2e20.10e3)') C_MolType(Tm),
& (-1.0d0/beta + H_test_1(Ib,Tm)/H_test_2(Ib,Tm) - H_test_3(Ib,Tm)/Counter, Ib=1,N_Box)
END DO
WRITE(6,*)
WRITE(6,'(A)') "Partial volume"
DO Tm=1,N_MolType
WRITE(6,'(1x,A19,18x,2e20.10e3)') C_MolType(Tm),
& (V_test_1(Ib,Tm)/H_test_2(Ib,Tm) - V_test_3(Ib,Tm)/Counter, Ib=1,N_Box)
END DO
WRITE(6,*)
WRITE(6,'(A)') "Heat Capacity"
DO Tm=1,N_MolType
WRITE(6,'(1x,A19,18x,2e20.10e3)') C_MolType(Tm),
& (1.0d0/(beta*beta) + (H_test_1(Ib,Tm)/H_test_2(Ib,Tm))*(H_test_1(Ib,Tm)/H_test_2(Ib,Tm))
& - H2_test_1(Ib,Tm)/H_test_2(Ib,Tm) - (H_test_3(Ib,Tm)/Counter)*(H_test_3(Ib,Tm)/Counter)
& + H2_test_3(Ib,Tm)/Counter , Ib=1,N_Box)
END DO
END IF
RETURN
END
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment