search

crest-lab / crest

CREST - A program for the automated exploration of low-energy molecular chemical space.
https://crest-lab.github.io/crest-docs/
GNU Lesser General Public License v3.0
199 stars 42 forks source link

Issue with initial geometry when using -nci mode #98

Closed BradenDKelly closed 2 years ago

BradenDKelly commented 2 years ago

Hello,

I get the following error when trying to use -nci mode on pairs of dimers.

Here is the command line of the call and error:

(base) [bkelly@c110370 paracetamol_anthracene]$ sudo crest enmin.pdb -nci

       ==============================================
       |                                            |
       |                 C R E S T                  |
       |                                            |
       |  Conformer-Rotamer Ensemble Sampling Tool  |
       |          based on the GFN methods          |
       |             P.Pracht, S.Grimme             |
       |          Universitaet Bonn, MCTC           |
       ==============================================
       Version 2.11.2, Fr 17. Dec 12:10:44 CEST 2021
  Using the xTB program. Compatible with xTB version 6.4.0

   Cite work conducted with this code as

   P. Pracht, F. Bohle, S. Grimme, PCCP, 2020, 22, 7169-7192.

   and  S. Grimme, JCTC, 2019, 15, 2847-2862.

   with help from:
   C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
   C. Plett, P.Pracht, S. Spicher

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 Command line input:
 > crest enmin.pdb -nci

  -nci  : Special NCI mode for non-covalently bound complexes or clusters.
  Automatically generated ellipsoide potential for NCI mode:
> $wall
>   potential=polynomial
>   ellipsoid: 1851.763096819424, 1850.270189733211, 1838.511887320862, all

 -------------------------
 xTB Geometry Optimization
 -------------------------

  Initial geometry optimization failed!
  Please check your input.

Here is a sample pdb file:

TITLE     paracetamol_anthracene
REMARK    THIS IS A SIMULATION BOX
CRYST1 3000.000 3000.000 3000.000  90.00  90.00  90.00 P 1           1
MODEL        1
ATOM      1  C1  API     1     999.095 997.735 997.920  1.00  0.00
ATOM      2  C2  API     1    1000.187 997.442 998.929  1.00  0.00
ATOM      3  O3  API     1    1000.721 996.353 998.991  1.00  0.00
ATOM      4  N4  API     1    1000.514 998.502 999.746  1.00  0.00
ATOM      5  C5  API     1    1001.518 998.5491000.745  1.00  0.00
ATOM      6  C6  API     1    1002.249 999.7351000.916  1.00  0.00
ATOM      7  C7  API     1    1003.240 999.8221001.908  1.00  0.00
ATOM      8  C8  API     1    1003.509 998.7201002.729  1.00  0.00
ATOM      9  C9  API     1    1002.784 997.5311002.566  1.00  0.00
ATOM     10  C10 API     1    1001.786 997.4471001.578  1.00  0.00
ATOM     11  O11 API     1    1004.477 998.8091003.680  1.00  0.00
ATOM     12  H12 API     1     999.327 998.653 997.370  1.00  0.00
ATOM     13  H13 API     1     998.139 997.855 998.440  1.00  0.00
ATOM     14  H14 API     1     999.013 996.904 997.211  1.00  0.00
ATOM     15  H15 API     1    1000.018 999.366 999.580  1.00  0.00
ATOM     16  H16 API     1    1002.0551000.5881000.274  1.00  0.00
ATOM     17  H17 API     1    1003.8051000.7401002.038  1.00  0.00
ATOM     18  H18 API     1    1002.987 996.6701003.195  1.00  0.00
ATOM     19  H19 API     1    1001.227 996.5231001.457  1.00  0.00
ATOM     20  H20 API     1    1004.574 997.9361004.093  1.00  0.00
ATOM     21  C1  XXX     2    1001.5851003.3731000.725  1.00  0.00
ATOM     22  C2  XXX     2    1001.0681002.8411001.917  1.00  0.00
ATOM     23  C3  XXX     2     999.9341002.0131001.884  1.00  0.00
ATOM     24  C4  XXX     2     999.3131001.7151000.658  1.00  0.00
ATOM     25  C5  XXX     2     998.1701000.8981000.624  1.00  0.00
ATOM     26  C6  XXX     2     997.5401000.613 999.400  1.00  0.00
ATOM     27  C7  XXX     2     996.383 999.815 999.367  1.00  0.00
ATOM     28  C8  XXX     2     995.757 999.533 998.142  1.00  0.00
ATOM     29  C9  XXX     2     996.2851000.048 996.947  1.00  0.00
ATOM     30  C10 XXX     2     997.4371000.850 996.976  1.00  0.00
ATOM     31  C11 XXX     2     998.0671001.136 998.201  1.00  0.00
ATOM     32  C12 XXX     2     999.2161001.946 998.234  1.00  0.00
ATOM     33  C13 XXX     2     999.8351002.245 999.460  1.00  0.00
ATOM     34  C14 XXX     2    1000.9711003.074 999.497  1.00  0.00
ATOM     35  H15 XXX     2    1002.4621004.0141000.753  1.00  0.00
ATOM     36  H16 XXX     2    1001.5461003.0681002.867  1.00  0.00
ATOM     37  H17 XXX     2     999.5411001.6041002.811  1.00  0.00
ATOM     38  H18 XXX     2     997.7651000.4971001.549  1.00  0.00
ATOM     39  H19 XXX     2     995.967 999.4171000.288  1.00  0.00
ATOM     40  H20 XXX     2     994.861 998.918 998.118  1.00  0.00
ATOM     41  H21 XXX     2     995.798 999.832 996.000  1.00  0.00
ATOM     42  H22 XXX     2     997.8361001.252 996.049  1.00  0.00
ATOM     43  H23 XXX     2     999.6171002.356 997.311  1.00  0.00
ATOM     44  H24 XXX     2    1001.3761003.489 998.578  1.00  0.00
TER
ENDMDL

If I use instead an xyz file I get the following error:

(base) [bkelly@c110370 paracetamol_anthracene]$ sudo crest enmin.xyz -nci

       ==============================================
       |                                            |
       |                 C R E S T                  |
       |                                            |
       |  Conformer-Rotamer Ensemble Sampling Tool  |
       |          based on the GFN methods          |
       |             P.Pracht, S.Grimme             |
       |          Universitaet Bonn, MCTC           |
       ==============================================
       Version 2.11.2, Fr 17. Dec 12:10:44 CEST 2021
  Using the xTB program. Compatible with xTB version 6.4.0

   Cite work conducted with this code as

   P. Pracht, F. Bohle, S. Grimme, PCCP, 2020, 22, 7169-7192.

   and  S. Grimme, JCTC, 2019, 15, 2847-2862.

   with help from:
   C.Bannwarth, F.Bohle, S.Ehlert, S.Grimme,
   C. Plett, P.Pracht, S. Spicher

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 Command line input:
 > crest enmin.xyz -nci

  -nci  : Special NCI mode for non-covalently bound complexes or clusters.
  Automatically generated ellipsoide potential for NCI mode:
> $wall
>   potential=polynomial
>   ellipsoid: 25.36194351795649, 21.38243357306829, 19.93121506887154, all

 -------------------------
 xTB Geometry Optimization
 -------------------------
 Geometry successfully optimized.

------------------------------------------------
Generating MTD length from a flexibility measure
------------------------------------------------
 System flexiblity is set to 1.0 for NCI mode
 flexibility measure :   1.000
 t(MTD) / ps    :    24.0
 Σ(t(MTD)) / ps :   144.0 (6 MTDs)

-------------------------------------
Starting a trial MTD to test settings
-------------------------------------
 Trial MTD 1 did not converge!
 Reducing the time step to 1 fs and trying again...

 Trial MTD 2 did not converge!
 Time step can`t be reduced further!
 Trying to turn SHAKE off...

 Trial MTD 3 did not converge!
 Time step can`t be reduced further!
 Trying to turn SHAKE off...

 Trial MTD 4 did not converge!
 Time step can`t be reduced further!
 Trying to turn SHAKE off...

 Trial MTD 5 did not converge!
 Time step can`t be reduced further!
 Trying to turn SHAKE off...

 Trial MTD 6 did not converge!
 Automatic xtb restart failed 6 times!
 Please try other settings by hand.

And here is the xyz file:

44
paracetamol_anthracene
C   999.0950000000  997.7350000000  997.9200000000
C  1000.1870000000  997.4420000000  998.9290000000
O  1000.7210000000  996.3530000000  998.9910000000
N  1000.5140000000  998.5020000000  999.7460000000
C  1001.5180000000  998.5490000000 1000.7450000000
C  1002.2490000000  999.7350000000 1000.9160000000
C  1003.2400000000  999.8220000000 1001.9080000000
C  1003.5090000000  998.7200000000 1002.7290000000
C  1002.7840000000  997.5310000000 1002.5660000000
C  1001.7860000000  997.4470000000 1001.5780000000
O  1004.4770000000  998.8090000000 1003.6800000000
H   999.3270000000  998.6530000000  997.3700000000
H   998.1390000000  997.8550000000  998.4400000000
H   999.0130000000  996.9040000000  997.2110000000
H  1000.0180000000  999.3660000000  999.5800000000
H  1002.0550000000 1000.5880000000 1000.2740000000
H  1003.8050000000 1000.7400000000 1002.0380000000
H  1002.9870000000  996.6700000000 1003.1950000000
H  1001.2270000000  996.5230000000 1001.4570000000
H  1004.5740000000  997.9360000000 1004.0930000000
C  1001.5850000000 1003.3730000000 1000.7250000000
C  1001.0680000000 1002.8410000000 1001.9170000000
C   999.9340000000 1002.0130000000 1001.8840000000
C   999.3130000000 1001.7150000000 1000.6580000000
C   998.1700000000 1000.8980000000 1000.6240000000
C   997.5400000000 1000.6130000000  999.4000000000
C   996.3830000000  999.8150000000  999.3670000000
C   995.7570000000  999.5330000000  998.1420000000
C   996.2850000000 1000.0480000000  996.9470000000
C   997.4370000000 1000.8500000000  996.9760000000
C   998.0670000000 1001.1360000000  998.2010000000
C   999.2160000000 1001.9460000000  998.2340000000
C   999.8350000000 1002.2450000000  999.4600000000
C  1000.9710000000 1003.0740000000  999.4970000000
H  1002.4620000000 1004.0140000000 1000.7530000000
H  1001.5460000000 1003.0680000000 1002.8670000000
H   999.5410000000 1001.6040000000 1002.8110000000
H   997.7650000000 1000.4970000000 1001.5490000000
H   995.9670000000  999.4170000000 1000.2880000000
H   994.8610000000  998.9180000000  998.1180000000
H   995.7980000000  999.8320000000  996.0000000000
H   997.8360000000 1001.2520000000  996.0490000000
H   999.6170000000 1002.3560000000  997.3110000000
H  1001.3760000000 1003.4890000000  998.5780000000

What I have found does work, is first running crest enmin.xyz -omd, and follow this with crest crest_best.xyz -nci, but I would think I should be able to just skip to the nci. It makes me feel like something else is going on that I should make sure is in order.

For further information: the pdb files are generated by gromacs, that is what makes and does a basic steepest descent energy minimization on the dimer pairs, and it then converts the .g96 file to a pdb (and does a bad job of it, I have a script to add in the 77-78 column the element names when I then use iodata to convert the pdb to xyz.

Finally, the coordinates are so large because gromacs 2021 won't run proper vacuum, so I put the dimers in a really large box and run it with PBC (the images are so far away, they will never be interacted with). Perhaps it is the large coordinates throwing off the -nci call? but then, I would expect that to also mess up the -omd call too.

BradenDKelly commented 2 years ago

I just noticed that the two $wall's is completely different for the xyz and pdb file, despite them being the same coordinates

drmewes commented 2 years ago

Dear Braden,

I am pretty sure the issue is related to your starting structure. Have you tried optimizing the starting structure with GFN2-xTB first? Then you can entirely avoid GROMACs and coordinate conversion as a source of error. Also, sometimes the autogenerated $wall is too small, which leads to unstable simulations (like in your second example). Here it helps to use -wscal 1.2 or 1.4 to increase the size of the cavity. Edit: Since you have plain organic molecules, I would at least run the tests using -gfnff as this will be ~1000 times faster (or at least -gfnff//gfn2). Cheers, Jan

awvwgk commented 2 years ago

Please, never start user-space programs with elevated privileges. crest is going to create a whole bunch of files and subprocesses using the same privilege level.

Note that wall potentials are always centered on the origin, for the NCI mode the center of mass of your molecule should be close to the actual coordinate system origin.

BradenDKelly commented 2 years ago

@drmewes @awvwgk I will center the molecules, I was curious about this, which is why I mentioned it. That could explain things. I don't plan on using gromacs for the initial starting coordinates in the future, but since i have several thousand... so it goes. I mostly just use it to make the a dimer that doesn't have substantial atomic overlap. My suspicion was that using crest enmin.xyz -omd was essentially just centering the dimer, with the added burden of doing a conformation search as well. i noticed the coordinates in crest_best.xyz were near the origin.

The use of sudo is necessary at the moment and I too find it unsavory.

pprcht commented 2 years ago

I'll second @awvwgk's comment. The problem here are the shifted coordinates of your system with the wall potential close to the origin. You can see that the calculations crash in the xtb.out of the trial MDs. The different wall potential sizes for pdb and xyz input formats could be a bug (probably the pdb coordinates are not shifted to the origin before determination of the wall potential). For the xyz input the generated wall potential size looks about right, but if you observe any convergence problems refer to @drmewes suggestion of the -wscal command first.

pprcht commented 2 years ago

Also, there is no -omd command in crest.

BradenDKelly commented 2 years ago

Also, there is no -omd command in crest.

This is true, and kind of funny. I lose track of when I am using CREST or xTB. Telling CREST to do -omd just triggers xTB. CREST is smarter than me, in that regard

It does generate errors though, but after awhile it gets around to using xTB to do -omd

BradenDKelly commented 2 years ago

Thanks fellows, it works fine now when I center at 0 first. It is interesting that gfn2 predicts pi-pi stacking as the best dimer and gfnff predicts stacked, but then paracetamol rotated 90 degrees so its edge is against the flat face of anthracene. Which is best has to be answered using DFT or CCSD(T).

Any plans on making CENSO work with freeware QM programs, like PSI4?

Cheers,

Braden

awvwgk commented 2 years ago

Any plans on making CENSO work with freeware QM programs, like PSI4?

Orca is free for academic use and supported in CENSO.

BradenDKelly commented 2 years ago

Any plans on making CENSO work with freeware QM programs, like PSI4?

Orca is free for academic use and supported in CENSO.

I am not academic :(

awvwgk commented 2 years ago

Well, can't help with that.

Of course, you are welcome to contribute Psi4 support to CENSO.

BradenDKelly commented 2 years ago

Well, can't help with that.

Of course, you are welcome to contribute Psi4 support to CENSO.

I may be able to pull that off. It has been awhile since I got SIGSEGV.