Messages

First of all you have to look at the number of adsorbate molecules. A monte carlo cycle in raspa is defined as min(20, N) where N is the number of adsorbate molecules. So if there are more molecules in the simulation cell than you simply have more moves per cycle.

Then you say, that computational expense of one cycle should scale approx. with N^2 because you have to evaluate the energies.

Does this help?

I just solved my problem.

When entering a fugacity coefficient the Molfraction still has to be specified otherwise raspa uses default values and sets them to equimolar fractions. I thougt that the molefraction is already represented in the fugacity coefficient but since the reservoir is an ideal gas this should not be the case. Anyway keep in mind to specify Molfractions when simulated mixture in GCMC simulation is not equimolar.

best Christopher

Dear RASPA fellows,

I want to do a Simulation in the Grand canonical ensemble in a COF. So the thermodynamic state is already chosen and I want to specify a fugacity coefficient that I calculated with PC-SAFT EOS.

The coefficient is already calculated in a mixture let's say 0.8 Molfrac methane and 0.2 Molfrac butane.

How do i correctly put in all necessary things. I think need:
ExternalPressure
FugacityCoefficient
MolFraction
IdealGasRosenbluthweight (for butane while methane is 1.)

Why do I have to specify the molfrac? Is this input fine now?
Best Christopher

Please take care and test if it really works. It was just an idea - I am not sure ;)

Just a strange idea:

would it work if you place the point charges that differ as an extra atom without charges with a fixed bond of 0.0 directly on top of the vdW bead?

Maybe you can try this.

Did you try to remove the line:
# pair truncated/shifted tailcorrections

and the empty line below? I would also remove the # at the end of the file.

Does it already fix it?

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!

Hi again,

oh sorry - seems like i mixed something up - default values are

CutOff VDW : 12.000000 (144.000000)
CutOff VDW switching on: 10.800000 (116.640000)
CutOff charge-charge : 12.000000 (144.000000)
CutOff charge-charge switching on: 7.800000 (60.840000)
CutOff charge-bonddipole : 12.000000 (144.000000)
CutOff charge-bondipole switching on: 8.400000 (70.560000)
CutOff bonddipole-bonddipole : 12.000000 (144.000000)
CutOff bonddipole-bondipole switching on: 9.000000 (81.000000)

This is from version 2.0.38

Hi!

in the raspa.pdf you can find some information about cutoffs on page 18. General information about the cutoff radius can be found in any molecular simulation book e.g. Frenkel and Smit. Cutoff means that the potential is only considered up to a certain value of radius.

In raspa for every cutoff value there are defaults, but you should really take care of the cutoff radius yourself because they should be considered as a part of your forcefield. If you use a different cutoff than you will probably receive different results!

You can get default values from the output file of your simulation (if you don't set any value for cutoffs). You get:
CutOff VDW : 14.000000 (196.000000)
CutOff VDW switching on: 12.600000 (158.760000)
CutOff charge-charge : 14.000000 (196.000000)
CutOff charge-charge switching on: 9.100000 (82.810000)
CutOff charge-bonddipole : 12.000000 (144.000000)
CutOff charge-bondipole switching on: 8.400000 (70.560000)
CutOff bonddipole-bonddipole : 12.000000 (144.000000)
CutOff bonddipole-bondipole switching on: 9.000000 (81.000000)

It's always a good idea to check your input with the output file.

Best 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

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

I forgot to mention, that I added

Intra14VDWScalingValue 0.0
Intra14ChargeChargeScalingValue 0.5

to the simulation.input file.

Best 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

RASPA determines the 1-4 neighbours by looking at the bending defined in the molecule file. Did you define all bendings?

Hello and thanks a bunch for your help!

It's really usefull and saves me a lot of extra time! Great!