periodic boundary conditions in MBX

[Niels Mikkelsen]

Hello developers!

I have been using the MB-nrg PEF for CO2/H2O mixtures and after my previous question about dipole moments everything has been working. The PEF has been very nice to work with so far, and has helped me quite a lot! I have only been working with structures that did not have PBC, but I now need to use PBC and I am experiencing some problems.

As a sanity check to see if PBC was working, I took an initial CO2/H2O structure that I had used for MD runs at 5 kelvin without PBC, and placed it in a box much larger than the structure and enabled PBC. That way, because of the large size of the box, the structure should behave the same way during a MD run at 5 kelvin even with PBC enabled, since the cell is so large that structure won't feel the mirrored structures. However, when running the MD, the structure immediately breaks apart, signaling to me that PBC does not work the way I think it does.

My settings are the following:

    mbx_config = {

    'box': [40.0, 0.0, 0.0, 0.0, 40.0, 0.0, 0.0, 0.0, 40.0],

    'twobody_cutoff': 100,

    'threebody_cutoff': 5,

    'max_n_eval_1b': 500,

    'max_n_eval_2b': 500,

    'max_n_eval_3b': 500,

    'dipole_tolerance': 1e-16,

    'dipole_max_it': 100,

    'dipole_method': 'cg',

    'alpha_ewald_elec': 0,

    'grid_density_elec': 2.5,

    'spline_order_elec': 6,

    'alpha_ewald_disp': 0,

    'grid_density_disp': 2.5,

    'spline_order_disp': 6,

    'ignore_2b_poly': [],

    'ignore_3b_poly': [],

    'ttm_pairs': []

    }

Here is a gif with an example of what I mean when I say the structure falls apart https://imgur.com/a/ZdIRTzM (apologies for the bad quality)

I have done other MD runs with PBC with a much smaller cell that fits the size of my structures, and there everything behaved the way I expected, but because I was unsure how to check that the PBC was indeed working, the only test I could think of was the one I described above that did not really work.

Is there something I am missing regarding how PBC works with the PEF? Is there anywhere I can read about how PBC works, or perhaps some examples?

Thanks in advance

Best, Niels

[Francesco Paesani]

Hi Niels,

Thank you for your interest in MBX. PBC in MBX works as in any other MD software. What MD engine are you using?

I am wondering if you can post here your initial configuration. The two-body cutoff in you json file is definitely too large for PBC simulations but I want to make sure that that is the only problem with your simulations.

Thank you,

Francesco

[Niels Mikkelsen]

Forgot to add what MD engine I am using in my previous message. I am using langevin dynamics implemented in ASE: https://wiki.fysik.dtu.dk/ase/ase/md.html

[Niels Mikkelsen]

Hi Francesco

Thank you for the fast answer! I have not really changed my settings since I started using the force field, I simply copied them from one of the examples on the Github. But I can definitely see why the two-body cutoff might be too large. Would you suggest I change it to a lower value?

Here is my initial configuration in .xyz format. If you prefer it in another format, please let me know. I am used to work in ASE with .traj files, so I am just assuming the .xyz format is good enough to understand my initial configuration.

120
Lattice="40.0 0.0 0.0 0.0 40.0 0.0 0.0 0.0 40.0" Properties=species:S:1:pos:R:3 pbc="T T T"
O       21.78634558      17.76774183      20.21930352
H       21.46852577      18.37156397      20.91019710
H       22.02180646      18.35877521      19.48275330
O       20.96138987      19.50692144      22.21857371
H       21.26827482      19.18017261      23.07014751
H       19.98785714      19.54893487      22.28919412
O       22.39065547      19.46597033      18.11923416
H       22.49560555      20.39902491      18.36241317
H       21.67292366      19.42775664      17.47026678
O       18.22270974      19.69175320      22.26971459
H       17.73958597      18.87347831      22.53032813
H       18.09075177      20.32149901      22.98532720
O       19.54511307      21.82666470      20.19977817
H       20.37303282      22.00144983      20.69024787
H       19.60321877      20.90350072      19.93415377
O       19.02965897      23.34884062      18.05888362
H       19.16105644      22.83360121      18.88392545
H       18.08291013      23.57087121      18.03615016
O       16.73016748      17.54164517      22.75216511
H       17.05523668      16.63336349      22.81419333
H       16.14194879      17.56394083      21.98276117
O       19.17985786      21.55561230      16.05952601
H       18.55557665      21.84355621      15.38596695
H       19.18480519      22.26285753      16.74454671
O       21.49296747      24.51818118      18.45101965
H       21.72551453      25.31800688      17.97563896
H       20.56604448      24.32158886      18.23030486
O       18.72786044      17.07844660      16.85121269
H       19.40735697      17.70870236      16.56628896
H       19.18212592      16.26258275      17.07348618
O       20.49039472      19.12872543      16.06387055
H       21.09518462      19.06156313      15.32133042
H       20.03428285      19.98307512      15.96540971
O       15.27053402      18.09745610      20.35303383
H       15.57260573      19.03493955      20.33311632
H       14.31318852      18.13290937      20.40175366
O       21.95715583      22.00720316      21.50730445
H       21.80363065      21.10166316      21.82007920
H       22.48915211      21.93331726      20.70071772
O       22.88416507      22.19492770      18.84858646
H       23.82415694      22.30894374      18.68580203
H       22.46442900      23.05534968      18.66026231
O       16.84334712      20.60777736      17.53725081
H       17.70177597      20.67646513      17.09623889
H       16.58007473      19.66514212      17.47537203
O       16.26089358      23.32182580      17.72335265
H       16.24495577      22.35387837      17.66883485
H       15.64669742      23.53371785      18.42984533
O       16.22011265      17.96277045      17.62497439
H       17.03726723      17.46186596      17.48717108
H       15.89052197      17.74276854      18.50566243
O       16.41528491      20.52220428      20.30297489
H       16.77028549      20.66258843      19.41083410
H       17.17280454      20.38493190      20.89320863
O       22.99857190      16.64032798      22.41117383
H       22.78654906      16.91079455      21.50192405
H       23.59289344      15.89267651      22.31724179
O       17.65162802      14.88764812      23.13773089
H       17.68822876      14.57343014      24.04331245
H       17.48378945      14.12284978      22.58415276
C       18.84124573      18.05255587      19.58543955
O       19.16886409      19.08813783      19.17060920
O       18.49307561      17.02711669      20.00144810
C       17.85756677      23.19743934      21.99947011
O       18.33281084      22.61431745      22.88751081
O       17.34932096      23.78314891      21.13767811
C       14.70324555      22.48968123      21.26469084
O       14.37831570      22.82327329      20.20052267
O       14.99736949      22.17774445      22.34370922
C       25.05981184      19.12492837      20.53708948
O       25.60901263      18.22192156      21.01879947
O       24.52896896      20.03999897      20.05711969
C       18.35150491      17.89395413      25.34547584
O       18.39458544      16.73552323      25.44277994
O       18.30426108      19.05030682      25.26927395
C       24.63483315      20.10592449      16.27913402
O       25.01965181      20.94317166      16.98596417
O       24.26957334      19.28146126      15.54682705
C       24.30273245      18.67379112      23.85572154
O       23.84349041      19.52042953      23.20099970
O       24.78762205      17.86440004      24.52941455
C       24.01448131      25.00403545      16.48830855
O       23.56081950      26.01232538      16.84807133
O       24.46390336      23.99915353      16.12215019
C       15.61790102      21.38391940      14.86843299
O       16.60556285      21.91487042      14.55722006
O       14.62386312      20.86314869      15.15613304
C       16.80218686      17.73344199      14.72613694
O       17.19169859      18.82901274      14.78804465
O       16.41189100      16.64661171      14.63072081
C       21.03868458      22.70935350      24.31027609
O       20.98151157      21.56249454      24.48807130
O       21.08464627      23.85934738      24.15897413
C       20.13921372      15.79273739      22.14089127
O       20.35248600      14.80890100      21.56173467
O       19.91784882      16.77492316      22.72056958
C       21.88716662      22.68206166      15.27622394
O       22.28778298      21.59454495      15.21199263
O       21.49294126      23.77459167      15.32289546
C       17.49213163      14.33257346      19.23679160
O       17.58376631      14.65717545      18.12729154
O       17.38550814      13.98767462      20.34175936
C       23.78466435      17.14984703      17.50881691
O       24.75229969      17.49842252      18.04771945
O       22.83161683      16.77999416      16.95559796
C       21.48078998      24.76270285      21.34693529
O       20.32752107      24.63568066      21.31279323
O       22.63104276      24.92418927      21.38307311
C       25.34091068      25.03352884      20.16884619
O       25.01479917      24.48421131      19.19641240
O       25.68681148      25.59373824      21.12265599
C       24.45421543      22.67039368      22.83895210
O       24.98436917      22.62263968      21.80533645
O       23.95870708      22.71513475      23.88715299
C       21.66446433      15.09700262      19.00722373
O       20.60295813      15.17730493      18.54017639
O       22.72879766      14.97657581      19.45497077
C       21.71393729      16.75299457      24.91602989
O       21.77268682      15.59796633      24.98574765
O       21.64160719      17.91590626      24.90073319

Best, Niels

torsdag den 13. april 2023 kl. 21.44.02 UTC+2 skrev Francesco Paesani:

[Francesco Paesani]

Hi Niels,

Please use the settings below and let me know if they work. The key changes are for  'alpha_ewald_elec' and  'alpha_ewald_disp'. For simulations in PBC, these values cannot be zero, which is correct for simulations without PBC, since they determine the number of k-vectors in the PME. Another important change is 'twobody_cutoff' since you may want to have the 2B cutoff to be smaller than half of your box. I also changed the 'threebody_cutoff' just to be consistent with the 3B cutoff for MB-pol. Finally, in PBC, you can use a larger value for  'dipole_tolerance', which will speed up the simulation. I changed it to 1e-8, which is a good compromise between speed and accuracy.

All the best,

Francesco

    mbx_config = {
    'box': [40.0, 0.0, 0.0, 0.0, 40.0, 0.0, 0.0, 0.0, 40.0],

    'twobody_cutoff': 9,
    'threebody_cutoff': 4.5,

    'max_n_eval_1b': 500,
    'max_n_eval_2b': 500,
    'max_n_eval_3b': 500,

    'dipole_tolerance': 1e-8,

    'dipole_max_it': 100,
    'dipole_method': 'cg',

    'alpha_ewald_elec': 0.60,
    'grid_density_elec': 2.5,
    'spline_order_elec': 6,
    'alpha_ewald_disp': 0.60,

[Marc Riera]

Hello Niels,

As a word of caution, there is a primitive version of CO2 (3b) that will be used by default unless you set:

'ignore_3b_poly': [ ["co2","co2","co2"]]

That three-body co2 version is not completed, but it helps. However, we have not tested or published anything with it and we do not recommend using it (see https://github.com/paesanilab/MBX/blob/master/examples/PEFs/035_co2-co2_mb-nrg_2bnb/JSON_files/mbx_pbc.json) 

In the future, if you want to run any of the PEFs implemented in MBX, you can just look at the examples folder and combine the corresponding JSON files. For now, the JSON file that Francesco sent you + the ignore 3b modification will allow you to use the correct PEF for CO2. 

Please let us know if there is any other issue!

Marc

[Niels Mikkelsen]

Hi Francesco and Marc

Once again thank you so much for all the help. I tried the MD with the settings above, and it now works as expected.

[Francesco Paesani]

Hi Niels,

Excellent! We look forward to seeing the results of your CO2/H2O simulations with the MB-nrg PEFs.

All the best,

Francesco