Messages

[a test simulation shows that the number of adsorbates in the restart file (seven, see the snippet attached below) differs from that reported at the last "PrintEvery" interval in the output file (nine, see the snippet attached below)]

Dear Prof. Dubbeldam and RASPA community,

The "Restart and Crash-Recovery" section of the RASPA manual mentions that "The restart file is written at 'PrintEvery' intervals.". According to this statement, I expect that the positions, velocities and forces of the adsorbate atoms recorded in the restart file to match the results printed at the last "PrintEvery" interval in the output file. However, a test simulation shows that the number of adsorbates in the restart file (seven, see the snippet attached below) differs from that reported at the last "PrintEvery" interval in the output file (nine, see the snippet attached below), which makes me confused. Could you please clarify what is the relationship between the configuration recorded by the restart file and the results documented in the output file?

Here are my simulation.input and snippets of restart and output files.

simulation.input

Code Select Expand

SimulationType                MonteCarlo
NumberOfCycles                1000
NumberOfInitializationCycles  1000
PrintEvery                    200

ContinueAfterCrash no

ChargeMethod                  Ewald
Forcefield                    ExampleMOFsForceField
UseChargesFromCIFFile yes
CutOffVDW                     10
CutOffChargeCharge            10
RemoveAtomNumberCodeFromLabel yes

Framework             0
FrameworkName         Cu-BTC
UnitCells             1 1 1
HeliumVoidFraction    0.745895
ExternalTemperature   273
ExternalPressure    2500000.0


Component 0 MoleculeName              N2
            MoleculeDefinition        ExampleDefinitions
            TranslationProbability    1.0
            RotationProbability       1.0
            ReinsertionProbability    1.0
            SwapProbability           1.0
            CreateNumberOfMolecules   0
restart_Cu-BTC_1.1.1_273.000000_2.5e+06

Code Select Expand

Components: 1 (Adsorbates 7, Cations 0)
========================================================================
Component 0 (N2)
Fractional-molecule-id component 0: -1
Lambda-factors component 0:  0.000000 0.000000 0.000000
Number-of-biasing-factors component 0: 21
Biasing-factors component 0:  0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
Maximum-CF-Lambda-change component 0: 0.500000
Maximum-CBCF-Lambda-change component 0: 0.500000

Maximum-translation-change component 0: 1.000000,0.487500,1.000000
Maximum-translation-in-plane-change component 0: 0.000000,0.000000,0.000000
Maximum-rotation-change component 0: 0.490874 1.472622 0.490874

Reactions: 0

Component: 0     Adsorbate    7 molecules of N2
------------------------------------------------------------------------
output_Cu-BTC_1.1.1_273.000000_2.5e+06.data

Code Select Expand

Current cycle: 800 out of 1000
========================================================================================================

Net charge: 1.0103e-13 (F: 1.0103e-13, A: 0, C: 0)
Current Box:  26.34300   0.00000   0.00000 [A]   Average Box:  26.34300   0.00000   0.00000 [A]
               0.00000  26.34300   0.00000 [A]                  0.00000  26.34300   0.00000 [A]
               0.00000   0.00000  26.34300 [A]                  0.00000   0.00000  26.34300 [A]
Box-lengths:   26.34300  26.34300  26.34300 [A] Average:  26.34300  26.34300  26.34300 [A]
Box-angles:   90.00000  90.00000  90.00000 [degrees] Average:  90.00000  90.00000  90.00000 [degrees]
Volume: 18280.82098 [A^3] Average Volume: 18280.82098 [A^3]

Loadings per component:
----------------------------------------------------------------------------------------------------------------------------------------------------
Component 0 (N2), current number of integer/fractional/reaction molecules: 9/0/0 (avg.   9.08489), density:  22.90147 (avg.  23.11749) [kg/m^3]
absolute adsorption:   9.00000 (avg.   9.08489) [mol/uc],   0.9302462147 (avg.   0.9390209050) [mol/kg],  26.0594337300 (avg.  26.3052433421) [mg/g]
                      20.8505160951 (avg.  21.0471917914) [cm^3 STP/g],   18.3237652830 (avg.  18.4966070141) [cm^3 STP/cm^3]
excess adsorption:    -0.1951236435 (avg.  -0.1102297608) [mol/uc],  -0.0201681145 (avg.  -0.0113934242) [mol/kg],  -0.5649790729 (avg.  -0.3191694608) [mg/g]
                      -0.4520476299 (avg.  -0.2553719336) [cm^3 STP/g],   -0.3972666493 (avg.  -0.2244249183) [cm^3 STP/cm^3]
----------------------------------------------------------------------------------------------------------------------------------------------------
Degrees of freedom: 45 0 45 0
Number of Framework-atoms:    624
Number of Adsorbates:           9 (9 integer, 0 fractional, 0 reaction)
Number of Cations:              0 (0 integer, 0 fractional, 0 reaction)

Current total potential energy:             -3908.4716378953 [K]  (avg.       -3915.5804048016)
Current Host-Host energy:                     0.0000000000 [K]  (avg.           0.0000000000)
Current Host-Adsorbate energy:            -3762.2229554409 [K]  (avg.       -3738.1974921185)
Current Host-Cation energy:                   0.0000000000 [K]  (avg.           0.0000000000)
Current Adsorbate-Adsorbate energy:        -146.2486824543 [K]  (avg.        -177.3829126831)
Current Cation-Cation energy:                 0.0000000000 [K]  (avg.           0.0000000000)
Current Adsorbate-Cation energy:              0.0000000000 [K]  (avg.           0.0000000000)



+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Finishing simulation
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

["FlexibleFramework"]

I do not think "FlexibleFramework" and "FrameworkDefinitions" have to be explicitly specified when simulating a rigid framework. See the RASPA manual examples of Adsorption isotherm of N2 in a metal-organic framework (MOF), Henry coefficients, enthalpy of adsorption and Adsorption isotherm of CO2 in Cu-BTC for a sample simulation.input file. In addition, you may also compare the input files of these two examples (Diffusion of benzene in rigid IRMOF-1 VS in flexible IRMOF-10) to see that "FlexibleFramework" and "FrameworkDefinitions" are not needed for rigid framework simulation.

In a sum, for your purpose of simulating a rigid structure, it is enough to write

Code Select Expand

Framework 0
FrameworkName calf-20
UnitCells 1.0 1.0 1.0

It would be best if you could upload your simulation files including force field definition (force_field.def or force_field_mixing_rules.def), pseudo_atoms.def, framework.cif, helium.def and simulation.input so that possible repliers can try to reproduce your error. Currently I would suggest checking the definition of He atom in the force_field.def, force_field_mixing_rules.def, pseudo_atoms.def and helium.def. RASPA documentation has provided an example about how to simulate the helium void-fraction of a structure. Maybe you can use that example as a reference.

['DensityProfile_methane.vtk']

I have tried your input file, and both the 'DensityProfile_methane.vtk' and 'COMDensityProfile_methane.vtk' output files can be successfully imported into iRASPA. I suggest checking the versions of RASPA and iRASPA that you have installed. The latest version of RASPA is v2.0.47, so I doubt whether you really installed RASPA v2.0.4 (released in 2016) or it is just a typo. The VTK files generated by such an old version of RASPA might not be readable by iRASPA. I am not sure of the exact version of RASPA I am using, but it was at least released after 2022. I downloaded iRASPA from the Windows store, which was released in 2018. By the way, the lines

Code Select Expand

Ensemble NVT
TimeStep 0.0005 in your input file seem unusual. From my understanding, these options are used for molecular dynamics simulations and should not be used in a Grand Canonical Monte Carlo simulation, so I commented these two lines when testing your input file.

[DensityProfile.vtk]

Are you sure your DensityProfile.vtk file is correct? It seems the file should contain totally 3375000 data points, but you just list 10 points. (Maybe you omit other data points for clarity?) Check /RASPA/RASPA2/examples/Visualization/FullCell/Density/MFI/hexane/RASPA/simulation.input for a sample input file. By the way, iRASPA can import the DensityProfile.vtk file generated by RASPA (with a newest version) to visualize the density plot, but I find it lacks some functions such as displaying the colorbar. Perhaps you can try iRASPA to see if it can generate the correct results.

In the gallery, then primitives, then IRMOF-1, you can see such an example.

Thank you for your reply. I do find the example. However, I am struggling to reproduce such a picture. The illustrations in the gallery are very compelling, but I think it may be difficult for a beginner like me to figure out how to draw those pictures without a tutorial. I was wondering if it would be possible for you to write a tutorial about how to draw the pictures in the gallery in the future. I understand that creating tutorials requires time and effort, but I believe this tutorial would be invaluable to the community.

Dear Prof. Dubbeldam,

In some literatures the cavity of a MOF is indicated by a sphere, can such an image be created by iRASPA? If not, could you please suggest me a suitable way to create such an image?

Thanks for your help.

Dear Prof. Dubbeldam,

In the manual via the help menu of iRASPA, I found that the distance between two atoms coule be measured by "alt + left mouse clicking on atoms 1 and 2". However, I could not use this function on my Windows desktop. Since the manual seems to be written for Macs, I wonder whether the measuring distance function is available on Windows.

Thans for your help.

[Segmentation Fault]

Dear Prof. Dubbeldam and RASPA community,

I would like to report a problem that I encountered. When I used the 'charge-equilibration' method to compute the charges of some relatively large MOFs, the simulation crashed with a Segmentation Fault.

Take ZIF-20.cif as an example. This cif file is in the original 'structures' directory of RASPA (~/RASPA/simulations/share/raspa/structures/cif), whose cell length a, b and c are all equal to 45.4725 Å.

My simulation.input file is as

Code Select Expand


SimulationType                MonteCarlo
NumberOfCycles                0
NumberOfInitializationCycles  0
PrintEvery                    100
RestartFile                   no

Forcefield                    ExampleMOFsForceField
CutOff                        12.8

ChargeFromChargeEquilibration    yes
ChargeEquilibrationPeriodic      yes
ChargeEquilibrationEwald         yes
SymmetrizeFrameworkCharges       no

Framework             0
FrameworkName         ZIF-20
UnitCells             1 1 1
ExternalTemperature   298.0
ExternalPressure      0.0


The simulation terminated with a Segmentation Fault.

29555 Segmentation fault      (core dumped) ${RASPA_DIR}/bin/simulate $1

However, if I change the simulated MOF from ZIF-20 to Cu-BTC (simply modify the 'FrameworkName' in the above simulation.input file), whose cell length a, b and c are all equal to 26.343 Å, the simulation can run correctly.

In addition, if I run the simulation of Cu-BTC with 2*2*2 unit cells rather than 1*1*1 unit cell, the simulation will also crash with a Segmentation Fault.

A similar problem had been reported previously in this forum by neumannrf:

https://forums.iraspa.org/index.php?topic=1018.msg1252#msg1252


neumannrf simulated TER.cif (cell length a, b and c are 9.8070 Å,  23.6460 Å and 20.2420 Å, respectivly) with 3*2*2 unit cells and got a Segmentation Fault, but when he used 1*1*1 unit cell, the calculation ran correctly. This is similar to the case of Cu-BTC, but for ZIF-20, even running the simulation with 1*1*1 unit cell leads to a Segmentation Fault.

Following neumannrf's debug steps, I also recompiled RASPA (v2.0.47) with CFLAGS="-w -ggdb -O0" and executed (in the case of ZIF-20 with 1*1*1 unit cell)

Code Select Expand


gdb ~/RASPA/simulations/bin/simulate

(gdb) run
Starting program: ~/RASPA/simulations/bin/simulate
_cell_length_a: 45.472500
_cell_length_b: 45.472500
_cell_length_c: 45.472500
_cell_length_alpha: 90.000000
_cell_length_beta: 90.000000
_cell_length_gamma: 90.000000
_symmetry_space_group_name_Hall: P 1 found space group: 1
_symmetry_space_group_name_H-M: P 1 found space group: 1
_symmetry_Int_Tables_number: 1
space group found from symmetry elements: 1 (nr elements: 1)
End reading cif-file
'force_field.def' file not found and therefore not used

Program received signal SIGSEGV, Segmentation fault.
0x00007ffffe82d2f4 in PrecomputeFixedEwaldContributions () at ewald.c:1078
1078          Eikx[j*MaxNumberOfCoulombicSites+i].re=Eikx[(j-1)*MaxNumberOfCoulombicSites+i].re*Eikx[MaxNumberOfCoulombicSites+i].re-

(gdb) bt
#0  0x00007ffffe82d2f4 in PrecomputeFixedEwaldContributions () at ewald.c:1078
#1  0x00007ffffeb2b015 in ReadInput (
    input=0x8009b90 "SimulationType", ' ' <repeats 16 times>, "MonteCarlo\nNumberOfCycles", ' ' <repeats 16 times>, "0\nNumberOfInitializationCycles  0\nPrintEvery", ' ' <repeats 20 times>, "100\nRestartFile", ' ' <repeats 19 times>, "no\n\nForcefield", ' ' <repeats 17 times>...) at input.c:9030
#2  0x00007ffffeae35b0 in ReadInputFile (filename=0x8008910 "simulation.input") at input.c:208
#3  0x00007fffff117c7f in run (inputData=0x8008910 "simulation.input", inputCrystal=0x8008930 "",
    raspaDir=0x7ffffffee268 "~/RASPA/simulations/", stream=false) at run.c:97
#4  0x00000000080014ba in main (argc=1, argv=0x7ffffffedf28) at main.c:106

(gdb) frame 0
#0  0x00007ffffe82d2f4 in PrecomputeFixedEwaldContributions () at ewald.c:1078
1078          Eikx[j*MaxNumberOfCoulombicSites+i].re=Eikx[(j-1)*MaxNumberOfCoulombicSites+i].re*Eikx[MaxNumberOfCoulombicSites+i].re-

(gdb) p MaxNumberOfCoulombicSites
$1 = 1024
(gdb) p i
$2 = 4351
(gdb) p j
$3 = 9
(gdb) p  Eikx[j*MaxNumberOfCoulombicSites+i].re
Cannot access memory at address 0x7fff8c505000
(gdb) p Eikx[(j-1)*MaxNumberOfCoulombicSites+i].re
$4 = 0.99841316505325695


It seems the fault is due to the inaccessible memory, so I guess the charge equilibration method for relatively large MOFs might require a lot of memory. However, I observed that during the simulation, only abount 2 G memory was used whereas my computer has a total memory of 32 G. By the way, I ran the simulation in WSL (Windows Subsystem for Linux).

In a sum, my questions are: is the segmentation fault described above due to the implementation of the charge equilibration method in RASPA? Or is it because I do something wrong? How can I compute the charges of relatively large frameworks (such as ZIF-20) in RASPA?

Any suggestions would be greatly appreciated.

Best regards,
Lingxiao Qin