1_2_ethanediol TraPPE modell

[Christopher]

Dear Dr Dubbeldam,

I tried to do a GibbsEnsemble simulation of ethanediol with the TraPPE modell. First I had some problems reading in the IntraCoulomb interactions. In the Outputfile Raspa set a scaling value of e-314 for some of them. After setting some extra IntraVDW interactions that I just scaled to zero it now seems that raspa is reading them in properly. But the Simulation gets somehow stuck. There is no exchange between the boxes and also the gibbsvolumechangemove is not working anymore (is not accepted).  Here is my molecule.def file:

# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-]
0.0
0.0
0.0
# Number Of Atoms
6
# Number Of Groups
1
# Alkane-group
flexible
# number of atoms
6
# atomic positions
0 H_alc
1 O_alc
2 CH2_alc
3 CH2_alc
4 O_alc
5 H_alc
# Chiral centers Bond  BondDipoles Bend  UrayBradley InvBend  Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb
               0    7            0    4            0       0        3            0         0            0         0               0            0        6           6
# Bond stretch: atom n1-n2, type, parameters
0 1 FIXED_BOND 0.945
1 2 FIXED_BOND 1.43
2 3 FIXED_BOND 1.54
3 4 FIXED_BOND 1.43
4 5 FIXED_BOND 0.945
0 4 LJ_12_6_BOND 75000000 0
1 5 LJ_12_6_BOND 75000000 0
# Bond bending: atom n1-n2-n3, type, parameters
0 1 2 HARMONIC_BEND 55400 108.5
1 2 3 HARMONIC_BEND 50400 109.47
2 3 4 HARMONIC_BEND 50400 109.47
3 4 5 HARMONIC_BEND 55400 108.5
# Torsion n1-n2-n3-n4 type
0 1 2 3 TRAPPE_DIHEDRAL    0.0   209.82   -29.17    187.93
1 2 3 4 TRAPPE_DIHEDRAL  503.24    0.0   -251.62   1006.47
2 3 4 5 TRAPPE_DIHEDRAL    0.0   209.82   -29.17    187.93
# intra vdw
0 3
0 3
0 3
0 3
0 3
0 3
# intra C
0 3
0 4
0 5
1 4
1 5
2 5
# Number of config moves
0

I did exactly the same simulation for 1-5 pentanediol and everything is working fine (Same simulation.input, pseudo_atoms and forcefield files). This one does not have the LJ_12_6_BOND repulsive term. Can you point out anything that I missed?

Thanks in advance!
Best
Christopher

[Christopher]

I forgot to mention, that I added

Intra14VDWScalingValue 0.0
Intra14ChargeChargeScalingValue 0.5

to the simulation.input file.

Best Christopher

[David Dubbeldam]

Are you using the latest version? (2.0.37).
I see no immediate issues with it.

        number of intra VDW interactions: 6
        --------------------------------------------
        interaction 0: A=0 B=3  scaling factor: 0
        interaction 1: A=0 B=3  scaling factor: 0
        interaction 2: A=0 B=3  scaling factor: 0
        interaction 3: A=0 B=3  scaling factor: 0
        interaction 4: A=0 B=3  scaling factor: 0
        interaction 5: A=0 B=3  scaling factor: 0

        number of intra charge-charge Coulomb interactions: 6
        -----------------------------------------------------------
        interaction 0: A=0 B=3  scaling factor: 0.5
        interaction 1: A=0 B=4  scaling factor: 1
        interaction 2: A=0 B=5  scaling factor: 1
        interaction 3: A=1 B=4  scaling factor: 0.5
        interaction 4: A=1 B=5  scaling factor: 1
        interaction 5: A=2 B=5  scaling factor: 0.5

Are you getting this?

[Christopher]

Dear Dr. Dubbeldam:
thanks for taking the time to help me!
I'm getting the same output. And also the forcefield in the output file looks fine. The version in the outputfile is
RASPA 2.0.38
Compiled as a 64-bits application
Compiler: gcc 7.4.1 20190905 [gcc-7-branch revision 275407]
Compile Date = Jul  1 2020, Compile Time = 13:22:41

I use this version of force_field_mixing_rules.def:
# general rule for shifted vs truncated
truncated
# general rule tailcorrections
yes
# number of defined interactions
13
# type interaction
He             LENNARD_JONES    10.9     2.64
CH4_sp3        LENNARD_JONES   148.0     3.73
CH3_sp3        LENNARD_JONES    98.0     3.75
CH2_sp3        LENNARD_JONES    46.0     3.95
CH_sp3         LENNARD_JONES    10.0     4.68
C_sp3          LENNARD_JONES     0.5     6.4
O_co2          LENNARD_JONES    79.0     3.05
C_co2          LENNARD_JONES    27.0     2.80
N_n2           LENNARD_JONES    36.0     3.31
O_alc          LENNARD_JONES    93.0     3.02
H_alc          LENNARD_JONES     0.0     0.0
CH2_alc        LENNARD_JONES    46.0     3.95
dummy          LENNARD_JONES     0.0     0.0
# general mixing rule for Lennard-Jones
Lorentz-Berthelot


And pseudo_atoms.def
#number of pseudo atoms
13
#type      print   as    chem  oxidation   mass        charge   polarization B-factor radii  connectivity anisotropic anisotropic-type   tinker-type
CH4_sp3    yes     C     C     0           16.04246    0.0      0.0          1.0      1.00   0            0           absolute           0
CH3_sp3    yes     C     C     0           15.03452    0.0      0.0          1.0      1.00   0            0           absolute           0
CH2_sp3    yes     C     C     0           14.02658    0.0      0.0          1.0      1.00   0            0           absolute           0
CH_sp3     yes     C     C     0           13.01864    0.0      0.0          1.0      1.00   0            0           absolute           0
C_sp3      yes     C     C     0           12.0        0.0      0.0          1.0      1.00   0            0           absolute           0
O_alc      yes     O     O     0           15.9994    -0.7      0.0          1.0      1.00   0            0           relative           0
H_alc      yes     H     H     0            1.00794    0.435    0.0          1.0      1.00   0            0           relative           0
CH2_alc    yes     C     C     0           14.02658    0.265    0.0          1.0      1.00   0            0           relative           0
C_co2      yes     C     C     0           12.0107     0.7      1.508        1.0      0.720  0            0           relative           0
O_co2      yes     O     O     0           15.9994    -0.35     0.9475       1.0      0.68   0            0           relative           0
N_n2       yes     N     N     0           14.00674   -0.482    0.0          1.0      0.7    0            0           relative           0
N_com      no      N     N     0           0.0         0.964    0.0          1.0      0.7    0            0           relative           0
dummy      no      N     N     0           0.0         0.0      0.0          1.0      1.00   0            0           relative           0

Have you ever tried TraPPE models of diols? I still can't figure out any reason, why it is not working.
If you have any hint, I would really appreciate!

Best
Christopher

[Christopher]

Maybe simulation.input would be interessting:

SimulationType                    MonteCarlo
NumberOfCycles                    50000
NumberOfInitializationCycles      50000
PrintEvery                        10
RestartFile                       no

CutOff                            14.0
CutOffVDW                         14.0
CutOffCoulomb                     14.0
EwaldPrecision                    1e-12
Forcefield                        local
WriteBinaryRestartFileEvery       1000
ContinueAfterCrash                no

Box 0
BoxLengths 40 40 40
BoxAngles 90 90 90
ExternalTemperature 600
ComputeMolecularPressure yes

Box 1
BoxLengths 40 40 40
BoxAngles 90 90 90
ExternalTemperature 600
ComputeMolecularPressure yes

GibbsVolumeChangeProbability 0.1

Component 0 MoleculeName                      1_2_ethanediol
            StartingBead                      0
            MoleculeDefinition                local
            Intra14ChargeChargeScalingValue   0.5
            Intra14VDWScalingValue            0
            TranslationProbability            1.0
            RandomTranslationProbability      1.0
            RotationProbability               1.0
            GibbsSwapProbability              1.0
            ReinsertionProbability            1.0
            CBMCProbability                   1.0
            CreateNumberOfMolecules           100 100


I varied EwaldPrecision from 1e-6 to 1e-12 and also other starting conditions like box lenghts or the number of molecules in a box. But everytime the simulation get's stuck.

Best Christopher

[dubbelda]

I copied your input, and see no problems.
Have you tried to run it with zero cycles? Does it finish? is the output precisely (pseudo-atoms, force field etc) exactly what you put in?
If that works, run it with 10 cycles and see if it works.

One thing that could go wrong is that your input files are not ascci-files, but utf8 or windows-files that contain different newlines.
With utilities like 'dos2unix' and 'unix2dos' you convert between windows and linux/Mac.

[Christopher]

I checked your remarks. The Output is the same that I put in. The Simulation runs as well with 0 as with 10 cycles. File formats are fine as well. I'm working solely on linux systems and i double checked it.

I now did the following:
In the molecule file I removed:
0 4 LJ_12_6_BOND 750000000 0
1 5 LJ_12_6_BOND 750000000 0
Then the molecule looks fine in the Movies.

If i include this part the bondlengths which are supposed to be fixed are not fixed anymore. So I thougt there is a problem with energyoverlap and set
EnergyOverlapCriteria             10e300        (the Energys at the Bondlengths become something around e+70)
MinimumRosenBluthFactor      10e-300     or          MinimumRosenBluthFactor      0
Still molecules are not grown correctly.

Do you have any other idea? I would really appreciate!