Periodic RESP charges for paracetamol form I crystal

514 views

Skip to first unread message

mbellucci

unread,

Mar 1, 2014, 3:02:12 PM3/1/14

to cp...@googlegroups.com

Dear CP2K Users,

I am trying to calculate the RESP charges for a unit cell of the form I crystal of paracetamol using periodic boundary conditions using version 2.5 of CP2K. I have tried many different things but in every calculation, the RESP charges seem to be unphysical, i.e. the magnitudes of the charges seem to very large and most have the wrong sign. The things I have tried in order to fix this include:

increasing the cutoff

changing the functional

increasing the size of the basis set

changing the values of the LENGTH and RANGE keywords in the PERIODIC_SYS section of RESP

changing the SURF_DIRECTION in the PERIODIC_SYS section of RESP

I've also experimented with the keywords WIDTH, RESTRAIN_HEAVIES_STRENGTH, and RESTRAIN_HEAVIES_TO_ZERO in RESP

Nothing seems to work for the RESP charges. However, when I turn off the periodic boundary conditions the calculations produce RESP charges that make physical sense. Also, for calculations with or without periodic boundary conditions, the DDAP charges make physical sense, but I'm hesitant to use them in MD simulations since they are less frequently used in MD simulations in the literature. If anyone could help identify an error I am making when calculating these RESP charges, I would greatly appreciate it! Also, if anyone could comment on the reliability of DDAP charges in general, I would appreciate it as well.

For the calculation, the coordinates of the form I crystal as well as the dimensions of the unit cell were obtained from the cambridge structural database. The input file I'm using, the outputted RESP charges, and the outputted DDAP charges are the following:

&GLOBAL

PROJECT form_1

RUN_TYPE ENERGY_FORCE

PRINT_LEVEL MEDIUM

WALLTIME 86400

&END GLOBAL

&FORCE_EVAL

METHOD Quickstep

&DFT

BASIS_SET_FILE_NAME EMSL_BASIS_SETS

POTENTIAL_FILE_NAME POTENTIAL

CHARGE 0

MULTIPLICITY 1

&MGRID

CUTOFF 300

REL_CUTOFF 50

&END MGRID

&QS

METHOD GAPW

EPS_DEFAULT 1.0E-12

MAP_CONSISTENT TRUE

EXTRAPOLATION ASPC

EXTRAPOLATION_ORDER 3

&END QS

&SCF

MAX_SCF 500

SCF_GUESS ATOMIC

EPS_SCF 1.0E-7

&OT ON

MINIMIZER DIIS

PRECONDITIONER FULL_ALL

ENERGY_GAP 0.001

&END OT

&OUTER_SCF

EPS_SCF 1.0E-5

MAX_SCF 100

&END OUTER_SCF

&END SCF

&XC

&XC_FUNCTIONAL PBE

&END XC_FUNCTIONAL

&END XC

&END DFT

&SUBSYS

&CELL

ABC 12.93 9.40 7.10

ALPHA_BETA_GAMMA 90.00 115.90 90.00

PERIODIC XYZ

&END CELL

&COORD

C 9.510000 3.334000 5.417000

C 10.770000 3.348000 4.827000

C 11.771000 2.497000 5.288000

C 11.518000 1.636000 6.343000

C 10.266000 1.639000 6.953000

C 9.279000 2.483000 6.489000

C 8.231000 4.763000 3.830000

C 6.944000 5.542000 3.720000

H 10.930000 3.957000 4.145000

H 12.675000 2.557000 4.918000

H 10.101000 1.062000 7.690000

H 8.407000 2.482000 6.891000

H 12.348000 0.517000 7.562000

H 7.724000 4.239000 5.582000

H 6.614000 5.480000 2.836000

H 6.512000 5.574000 4.484000

H 7.181000 6.589000 3.692000

N 8.426000 4.170000 5.017000

O 12.536000 0.800000 6.745000

O 9.035000 4.678000 2.901000

C 6.784000 8.034000 0.970000

C 5.524000 8.048000 1.560000

C 4.523000 7.197000 1.099000

C 4.776000 6.336000 0.044000

C 6.028000 6.339000 -0.567000

C 7.015000 7.183000 -0.102000

C 8.063000 9.463000 2.557000

C 9.350000 10.242000 2.667000

H 5.364000 8.657000 2.242000

H 3.618000 7.257000 1.469000

H 6.193000 5.762000 -1.303000

H 7.887000 7.182000 -0.505000

H 3.945000 5.217000 -1.175000

H 8.570000 8.939000 0.805000

H 9.680000 10.180000 3.551000

H 9.782000 10.274000 1.903000

H 9.113000 11.289000 2.695000

N 7.867000 8.870000 1.370000

O 3.758000 5.500000 -0.358000

O 7.259000 9.378000 3.486000

C 0.319000 6.066000 0.970000

C -0.941000 6.052000 1.560000

C -1.942000 6.903000 1.099000

C -1.689000 7.764000 0.044000

C -0.437000 7.761000 -0.567000

C 0.550000 6.917000 -0.102000

C 1.598000 4.637000 2.557000

C 2.885000 3.858000 2.667000

H -1.101000 5.443000 2.242000

H -2.847000 6.843000 1.469000

H -0.272000 8.338000 -1.303000

H 1.422000 6.918000 -0.505000

H -2.520000 8.883000 -1.175000

H 2.105000 5.161000 0.805000

H 3.215000 3.920000 3.551000

H 3.317000 3.826000 1.903000

H 2.648000 2.811000 2.695000

N 1.402000 5.230000 1.370000

O -2.707000 8.600000 -0.358000

O 0.794000 4.722000 3.486000

C 3.045000 1.366000 5.417000

C 4.305000 1.352000 4.827000

C 5.306000 2.203000 5.288000

C 5.053000 3.064000 6.343000

C 3.801000 3.061000 6.953000

C 2.814000 2.217000 6.489000

C 1.766000 -0.063000 3.830000

C 0.479000 -0.842000 3.720000

H 4.465000 0.743000 4.145000

H 6.210000 2.143000 4.918000

H 3.636000 3.638000 7.690000

H 1.942000 2.218000 6.891000

H 5.883000 4.183000 7.562000

H 1.259000 0.461000 5.582000

H 0.149000 -0.780000 2.836000

H 0.047000 -0.874000 4.484000

H 0.716000 -1.889000 3.692000

N 1.961000 0.530000 5.017000

O 6.071000 3.900000 6.745000

O 2.570000 0.022000 2.901000

&END COORD

&KIND H

BASIS_SET 6-31++Gxx

POTENTIAL GTH-PBE-q1

&END KIND

&KIND C

BASIS_SET 6-31++Gxx

POTENTIAL GTH-PBE-q4

&END KIND

&KIND N

BASIS_SET 6-31++Gxx

POTENTIAL GTH-PBE-q5

&END KIND

&KIND O

BASIS_SET 6-31++Gxx

POTENTIAL GTH-PBE-q6

&END KIND

&END SUBSYS

&PROPERTIES

&RESP

&PERIODIC_SYS

ATOM_LIST 1..80

LENGTH 3.0

RANGE 3.0 5.0

SURF_DIRECTION Z

&END PERIODIC_SYS

&PRINT

&V_RESP_CUBE

&END V_RESP_CUBE

&RESP_CHARGES_TO_FILE

&END RESP_CHARGES_TO_FILE

&PROGRAM_RUN_INFO

&END

&COORD_FIT_POINTS

&END

&END PRINT

&END RESP

&FIT_CHARGE

&END FIT_CHARGE

&END PROPERTIES

&END FORCE_EVAL

RESP charges:

Type | Atom | Charge

RESP 1 C 0.070078

RESP 2 C 0.903292

RESP 3 C 1.113722

RESP 4 C -0.120125

RESP 5 C 1.058323

RESP 6 C 0.655164

RESP 7 C 0.200902

RESP 8 C 2.666416

RESP 9 H -1.222912

RESP 10 H -1.170830

RESP 11 H -0.983979

RESP 12 H -0.818745

RESP 13 H -1.214944

RESP 14 H -1.366249

RESP 15 H -1.224385

RESP 16 H -1.119378

RESP 17 H -1.137280

RESP 18 N 1.259042

RESP 19 O 1.713353

RESP 20 O 0.738109

RESP 21 C 0.041413

RESP 22 C 0.956097

RESP 23 C 1.040182

RESP 24 C -0.082712

RESP 25 C 1.051482

RESP 26 C 0.595410

RESP 27 C 0.202160

RESP 28 C 2.662283

RESP 29 H -1.166123

RESP 30 H -1.126516

RESP 31 H -0.984189

RESP 32 H -0.800838

RESP 33 H -1.188989

RESP 34 H -1.401436

RESP 35 H -1.215709

RESP 36 H -1.113423

RESP 37 H -1.128041

RESP 38 N 1.271056

RESP 39 O 1.679915

RESP 40 O 0.705805

RESP 41 C 0.116651

RESP 42 C 0.813652

RESP 43 C 1.108160

RESP 44 C -0.087035

RESP 45 C 1.066794

RESP 46 C 0.598212

RESP 47 C 0.201601

RESP 48 C 2.665972

RESP 49 H -1.193708

RESP 50 H -1.119339

RESP 51 H -0.987625

RESP 52 H -0.794927

RESP 53 H -1.193920

RESP 54 H -1.403616

RESP 55 H -1.227170

RESP 56 H -1.113105

RESP 57 H -1.132964

RESP 58 N 1.272201

RESP 59 O 1.669491

RESP 60 O 0.737908

RESP 61 C 0.017005

RESP 62 C 0.961044

RESP 63 C 1.064429

RESP 64 C -0.116028

RESP 65 C 1.044203

RESP 66 C 0.665757

RESP 67 C 0.203199

RESP 68 C 2.660164

RESP 69 H -1.159143

RESP 70 H -1.166104

RESP 71 H -0.983609

RESP 72 H -0.830205

RESP 73 H -1.208912

RESP 74 H -1.380887

RESP 75 H -1.210223

RESP 76 H -1.114933

RESP 77 H -1.129467

RESP 78 N 1.264612

RESP 79 O 1.718065

RESP 80 O 0.706396

Total 0.000000

DDAP FULL DENSITY charges:

Atom | Charge

1 C 0.164768

2 C 0.038789

3 C -0.067304

4 C 0.359095

5 C -0.273509

6 C -0.132861

7 C 0.453357

8 C -0.073176

9 H -0.008341

10 H -0.007149

11 H 0.160840

12 H 0.127025

13 H 0.128967

14 H 0.094082

15 H -0.056650

16 H -0.027392

17 H -0.004180

18 N -0.249484

19 O -0.299466

20 O -0.327283

21 C 0.164803

22 C 0.039171

23 C -0.067954

24 C 0.359126

25 C -0.273851

26 C -0.132660

27 C 0.453316

28 C -0.072910

29 H -0.008512

30 H -0.006834

31 H 0.161121

32 H 0.126815

33 H 0.128910

34 H 0.094007

35 H -0.056723

36 H -0.027522

37 H -0.004288

38 N -0.249623

39 O -0.299372

40 O -0.327146

41 C 0.164786

42 C 0.039184

43 C -0.067967

44 C 0.359125

45 C -0.273855

46 C -0.132651

47 C 0.453282

48 C -0.072882

49 H -0.008520

50 H -0.006829

51 H 0.161118

52 H 0.126809

53 H 0.128911

54 H 0.094000

55 H -0.056734

56 H -0.027529

57 H -0.004294

58 N -0.249597

59 O -0.299364

60 O -0.327128

61 C 0.164787

62 C 0.038773

63 C -0.067297

64 C 0.359102

65 C -0.273505

66 C -0.132867

67 C 0.453392

68 C -0.073203

69 H -0.008330

70 H -0.007154

71 H 0.160838

72 H 0.127032

73 H 0.128964

74 H 0.094090

75 H -0.056640

76 H -0.027383

77 H -0.004176

78 N -0.249513

79 O -0.299471

80 O -0.327306

Total 0.000000

Dorothea Golze

unread,

Mar 1, 2014, 3:19:58 PM3/1/14

to cp...@googlegroups.com

Hi,

the keywords "RANGE" and "LENGTH" make no sense in your case. They are for slab-like periodic systems. The section "PERIODIC_SYS" and "NONPERIODIC_SYS" refer to the way the fitting points are collected.

PERIODIC_SYS sets the parameter for collecting fitting points for slab-like systems. In "NONPERIODC_SYS" fitting points are collected in spheres around the atoms. That is what you need here. Whether a periodic RESP fit is required or not is solely decided on the usage of a periodic or non-periodic Poisson solver.
I admit the names of these subsection are maybe a little bit misleading. I should change them....

There were questions about the periodic RESP charges before, see

my last post there:
https://groups.google.com/forum/#!msg/cp2k/G3THGbvvIcU/CgGZxvtbGrsJ

Kind regards,

Dorothea

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/groups/opt_out.

mbellucci

unread,

Mar 2, 2014, 2:35:53 PM3/2/14

to cp...@googlegroups.com, dorothe...@googlemail.com

Hi Dorothea,

Thanks for the reply. I tried what you suggested, i.e. using the NONPERIODIC_SYS keyword instead of the PERIODIC_SYS keyword, but it only seems to give physically meaningful charges when I make the cell dimensions very large (25.0 25.0 25.0 cubic box) so that all 4 molecules that make up the unit cell are completely within the box. However, these very large cubic cell dimensions don't correspond to the actual dimensions of the monoclinic unit cell, and so, the periodicity with this large cell doesn't correspond the periodicity of the actual crystal. When I try to use the cell dimensions of the unit cell (see previously posted input file) with the NONPERIODIC_SYS keyword, once again I end up with non-physical RESP charges. I've also tried a range of different values for RMIN and RMAX in the NONPERIODIC_SYS section, but this doesn't seem to help very much. Perhaps I am missing something simple, and so any suggestions would be appreciated.

Mike

Dorothea Golze

unread,

Mar 2, 2014, 4:35:55 PM3/2/14

to mbellucci, cp...@googlegroups.com

Hi Mike,

there are a few things you should check on

1) Search for the Root-mean-square (RMS) and Relative root-mean-square (RRMS) error in the output file. They give you hints on the quality of the fit. Check J . Phys. Chem. 1993,97, 10269-10280

to get an idea in which range the values should be for ESP.

2) You can directly compare the potentials (original and fitted one). Enable in the respective print sections &V_HARTREE_CUBE and &V_RESP_CUBE

Sometimes the fitting is good, but charges are unphysical. This can theoretically happen with ESP;

To make the charges "more physical" you have to employ restraints (restrained ESP -> RESP)

3) check if the fitting points were collected as you intended by switching on &COORD_FIT_POINTS in the print section. Fix the first line in the generated .xyz (number atoms + number fitting points) to view in vmd

4) if the collection of fitting points is correct and the fitting is bad according to the RRMS values, please send me your latest resp input file; It should work with monoclinic cells. However, I never explicitly tested it (only orthorhombic and hexagonal cells)

5) Sometimes the potential is "too complicated" and the fitting is just bad. This can also happen, even for small molecules, see J . Phys. Chem. 1993,97, 10269-10280

Kind regards,

Dorothea

hut...@chem.uzh.ch

unread,

Mar 3, 2014, 3:33:14 AM3/3/14

to cp...@googlegroups.com

Hi

please remember that ESP calculations work by fitting the
electrostatic field away from the charge distribution.
In a dense periodic system (closely packed crystal or a liquid)
there is no space without charge density and you cannot use
the normal method.
The periodic ESP code in CP2K is targeted for open systems
(zeolites, MOFs) and surfaces.
For liquids we developed a method based on Wannier functions
back in the CPMD days (J. Chem. Phys. 121, 5133-5142 (2004)),
not implemented in CP2K).
A work around would be to use QM/MM and self-consistently
get the charges.

regards

Juerg Hutter

--------------------------------------------------------------
Juerg Hutter                         Phone : ++41 44 635 4491
Institut für Chemie                  FAX   : ++41 44 635 6838
Universität Zürich                   E-mail: hut...@chem.uzh.ch
Winterthurerstrasse 190
CH-8057 Zürich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: mbellucci
Sent by: cp...@googlegroups.com
Date: 03/02/2014 08:36PM
Cc: dorothe...@googlemail.com
Subject: Re: [CP2K:5018] Periodic RESP charges for paracetamol form I crystal

mbellucci

unread,

Mar 3, 2014, 4:56:18 PM3/3/14

to cp...@googlegroups.com

Thank you Dorothea and Juerg for your helpful comments. I managed to get physically meaningful charges from the RESP fit by using restraints in the RESP fit and using the DDAP charges as the target values. At the moment I am varying the strength of the restraints in an effort to minimize the RMSD of the fit. Thank you for your help!

Best,

Mike

Reply all

Reply to author

Forward

0 new messages