DFT-D
Jörg Saßmannshausen
admittingly, my skills in cp2k are a bit rusty.
I wanted to look into the DFT-D option of cp2k, specially how the second
derivatives calculation are performing on rather large molecules.
From the manual I takt it one can do DFT-D calculations so I thought I am
using the water example and play around with that. I decided, mainly for
speed reasons, to go with the BP functional and compare that MD run with one
where DFT-D is switched on, i.e. BP-D.
However, I am not sure whether I done the setup correctly, so I thought I
might as well print out whether or not dispension is actually switch on.
Playing around with various option in
__ROOT__%FORCE_EVAL%DFT%PRINT%DFT_CONTROL_PARAMETERS
where somehow unsuccesful in this respect. Would somebody just briefly comment
on whether my input file (attached below) is correct and how to print out the
DFT parameters?
I have omitted the coordination section to save some space.
All the best from London!
Jörg
Input file:
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../BASIS_SET
POTENTIAL_FILE_NAME ../POTENTIAL
&MGRID
CUTOFF 280
&END MGRID
&QS
EPS_DEFAULT 1.0E-12
WF_INTERPOLATION PS
EXTRAPOLATION_ORDER 3
&END QS
&SCF
SCF_GUESS ATOMIC
&OT ON
MINIMIZER DIIS
&END OT
# SCF_GUESS RESTART
# EPS_SCF 1.0E-7
&PRINT
&RESTART OFF
&END
&END
&END SCF
&XC
&XC_FUNCTIONAL BP
&END XC_FUNCTIONAL
&VDW_POTENTIAL
&PAIR_POTENTIAL
Type Grimme
&end PAIR_POTENTIAL
# POTENTIAL_TYPE DISPERSION_FUNCTIONAL
&END VDW_POTENTIAL
&END XC
&print
&DFT_CONTROL_PARAMETERS high
&end DFT_CONTROL_PARAMETERS
&end print
&END DFT
&SUBSYS
&CELL
ABC 9.8528 9.8528 9.8528
&END CELL
# 32 H2O (TIP5P,1bar,300K) a = 9.8528
&COORD
[ ... ]
&END COORD
&KIND H
BASIS_SET TZV2P-GTH
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET TZV2P-GTH
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O-32-vdw
RUN_TYPE MD
PRINT_LEVEL low
&TIMINGS
THRESHOLD 0.000001
&END
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 300
TIMESTEP 0.5
TEMPERATURE 300.0
&END MD
&END MOTION
--
*************************************************************
Jörg Saßmannshausen
University College London
Department of Chemistry
Gordon Street
London
WC1H 0AJ
email: j.sassma...@ucl.ac.uk
web: http://sassy.formativ.net
Please avoid sending me Word or PowerPoint attachments.
See http://www.gnu.org/philosophy/no-word-attachments.html
hut...@pci.uzh.ch
here is an example for the Grimme D2 method
&vdW_POTENTIAL
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD2
R_CUTOFF 15.0
SCALING 1.0
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
The SCALING parameter refers to the s6 term in this method.
Default values for some functionals are available in the
code through
REFERENCE_FUNCTIONAL BLYP
This is an example for the new D3 method
DISPERSION_FUNCTIONAL PAIR_POTENTIAL
&PAIR_POTENTIAL
TYPE DFTD3
REFERENCE_FUNCTIONAL BLYP
CALCULATE_C9_TERM .TRUE.
PARAMETER_FILE_NAME dftd3.dat
R_CUTOFF 15.0
&END PAIR_POTENTIAL
&END vdW_POTENTIAL
With VERBOSE_OUTPUT TRUE you can get all kinds of detailed
energy contributions.
regards
Juerg Hutter
--------------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
---------------------------------------------------------------
-----cp...@googlegroups.com wrote: -----
To: cp2k <cp...@googlegroups.com>
From: Jörg Saßmannshausen <j.sassma...@ucl.ac.uk>
Sent by: cp...@googlegroups.com
Date: 11/21/2010 11:23AM
Subject: [CP2K:2928] DFT-D
Dear all
Jörg
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To post to this group, send email to cp...@googlegroups.com.
To unsubscribe from this group, send email to cp2k+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/cp2k?hl=en.
Matt W
> This is an example for the new D3 method
> DISPERSION_FUNCTIONAL PAIR_POTENTIAL
> &PAIR_POTENTIAL
> TYPE DFTD3
> REFERENCE_FUNCTIONAL BLYP
> CALCULATE_C9_TERM .TRUE.
> PARAMETER_FILE_NAME dftd3.dat
> R_CUTOFF 15.0
> &END PAIR_POTENTIAL
> &END vdW_POTENTIAL
>
Do you recommend the three body (C9) term? It seems Grimme was
somewhat dubious about it in the original JCP 2010 paper?
Cheers,
Matt
hut...@pci.uzh.ch
we do have analytic gradients for the C9 terms. Grimme was
negative in part because of the costs of their numerical
gradients. For small systems the C9 terms are certainly not important.
We do have the possibility to do C9 using fixed C6. This speeds
up the calculation also.
The following example shows the improvements of D3 for surface
systems, although I don't know if this includes the C9
http://www.scm.com/News/DFT-D3.html
With VERBOSE_OUTPUT you can check easily how big the contribution
of E9 is.
regards
Juerg
--------------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
---------------------------------------------------------------
-----cp...@googlegroups.com wrote: -----
To: cp2k <cp...@googlegroups.com>
From: Matt W <mattwa...@gmail.com>
Sent by: cp...@googlegroups.com
Date: 11/22/2010 02:42PM
Subject: [CP2K:2930] Re: DFT-D
Hi,
Cheers,
Matt
--
Jörg Saßmannshausen
many thanks for the prompt reply and the usefull informations.
I think I have to update my cp2k version which is from last summer ;-)
All the best from London
Jörg
--
Jörg Saßmannshausen
ok, I am now playing around with the DFT-D in cp2k.
I am currently trying to ascertain how much memory we would require for the
work we are planning to do and for the new cluster. Hence, these are test
runs of the right size of molecule but not with a really good optimised
structure (i.e. not optimised with DFT-D).
I have tried both the DFT-D2 and DFT-D3 option and I got 2 problems:
a) the test molecule needs about > 48 GB of RAM, which is significant. Is
there any way to get that into a more handable size? I have tried the
SAVE_MEM in GLOBAL but that did not really had any effect it seems. Following
Matthias' suggestion some time ago, I am using a cutoff of 660, which was
working well for different molecules in the past.
b) I get the following error for the DFT-D2 frequency calculation:
[ ... ]
REPLICA| layout of the replica grid, number of groups
2
REPLICA| layout of the replica grid, size of each group
1
REPLICA| MPI process to grid (group,rank) correspondence:
( 0 : 0, 0) ( 1 : 1, 0)
VIB| Vibrational Analysis Info
Pair potential vdW calculation
Dispersion potential type: DFTD2
Scaling parameter (s6) 1.00000000000000000
Exponential prefactor 20.000000000000000
Total vdW energy [au] : -6.37708902660472643E-002
Total vdW energy [kcal] : -40.016837470249470
Dispersion Forces
Atom Kind Forces
[ ... ]
|G| = 5.08859983571732719E-002
Stress Tensor (dispersion)
-0.911187517690E-03 -0.491992126627E-02 0.140501262854E-01
-0.491992126627E-02 0.134365275526E-01 0.179254944339E-02
0.140501262854E-01 0.179254944339E-02 -0.155716224075E-01
Tr(P)/3 : -1.01542745750825042E-003
CP2K| condition FAILED at line 253
CP2K| Abnormal program termination, stopped by process number 1
CP2K| condition FAILED at line 253
CP2K| Abnormal program termination, stopped by process number 0
For the DFT-D3 calculation I get this error:
[ ... ]
|G| = 1.19763995420629592E-002
Stress Tensor (dispersion)
-0.525748742588E-02 -0.260638877285E-02 0.799232865409E-02
-0.260638877285E-02 -0.433023621108E-03 0.181088946886E-02
0.799232865409E-02 0.181088946886E-02 -0.929520614835E-02
Tr(P)/3 : -4.99523906511300034E-003
CP2K| condition FAILED at line 236
CP2K| Abnormal program termination, stopped by process number 0
It appears to me I am doing something wrong. I basically have copied and
pasted your section into the force_eval.inc file.
Can you or somebody be so kind and point me in the right direction? It is
possible to do the frequency calculation with DFT-D or am I wrong here?
I have attached the input files dft-d-D2.inp, force_eval-D2.inc and
subsys.inc.
All the best from London
Jörg
On Montag 22 November 2010 hut...@pci.uzh.ch wrote:
Matthias Krack
your simulation box size seems to be about 27 Angstrom^3 and using a
cutoff 660 Ry for such a box will result in large grids and thus huge
memory allocations. Possibly, the extreme memory request is the reason
for abnormal program termination. You may check your grid sizes using
http://cp2k.berlios.de/manual/CP2K_INPUT/FORCE_EVAL/PRINT/GRID_INFORMATION.html
to get an idea of what will be used.
Afaik, SAVE_MEM helps only for systems which are large due to the
number of atoms employed. This is not the case for your system and
thus there is no impact from this keyword.
cheers,
Matthias
On Dec 13, 11:54 pm, Jörg Saßmannshausen <j.sassmannshau...@ucl.ac.uk>
wrote:
> > Seehttp://www.gnu.org/philosophy/no-word-attachments.html
> > --
> > You received this message because you are subscribed to the Google Groups
> > "cp2k" group. To post to this group, send email to cp...@googlegroups.com.
> > To unsubscribe from this group, send email to
> > cp2k+uns...@googlegroups.com. For more options, visit this group at
> >http://groups.google.com/group/cp2k?hl=en.
>
> --
> *************************************************************
> Jörg Saßmannshausen
> University College London
> Department of Chemistry
> Gordon Street
> London
> WC1H 0AJ
>
>
> Please avoid sending me Word or PowerPoint attachments.
> Seehttp://www.gnu.org/philosophy/no-word-attachments.html
>
> < 1KViewDownload
>
> subsys.inc
> 4KViewDownload
>
> force_eval-D2.inc
> 1KViewDownload
Jörg Saßmannshausen
thanks for the reply.
I was using the script from the workshop to get the box size. You are not
aware of any bugs in that script which could lead to a wrong box size for
rathe large but flat molecules?
I will toy around with the grid size and see if I can get it going.
Probably moving to a PBC calculation would not help with the problem at all I
would guess, or?
All the best from a rainy London
Jörg
email: j.sassma...@ucl.ac.uk
Matthias Krack
I don't know to which script you are referring. Anyhow, just take care
that there is at least a 3 Angstrom (better more, ie. 5 or 6 A)
distance between the atoms of your system and the box walls for
isolated system runs, e.g. the difference between the maximum and
minimum atomic x coordinate of your system plus 5 A will give you the
box edge length a (do the same for y and z to get b and c). The extra
space (vacuum portion) is needed, because decoupling methods like MT
require that the electron density is (almost) zero at the box walls.
Thus a flat molecule can be placed in an adapted flat box. Use a
visualization tool to check your setup. Note, for larger systems and
thus larger boxes (like 27 A**3) you will need proper computational
resources, i.e. a parallel computer with enough CPU cores and RAM to
keep the grids in memory. For your system you might need 32 CPU cores
or even more.
cheers,
Matthias
On Dec 16, 10:54 pm, Jörg Saßmannshausen <j.sassmannshau...@ucl.ac.uk>
wrote:
>
> thanks for the reply.
> I was using the script from the workshop to get the box size. You are not
> aware of any bugs in that script which could lead to a wrong box size for
> rathe large but flat molecules?
>
> I will toy around with the grid size and see if I can get it going.
> Probably moving to a PBC calculation would not help with the problem at all I
> would guess, or?
>
> All the best from a rainy London
>
> Jörg
>
> On Dienstag 14 Dezember 2010 Matthias Krack wrote:
>
> > Hi Joerg,
>
> > your simulation box size seems to be about 27 Angstrom^3 and using a
> > cutoff 660 Ry for such a box will result in large grids and thus huge
> > memory allocations. Possibly, the extreme memory request is the reason
> > for abnormal program termination. You may check your grid sizes using
>
> read more »
Teodoro Laino
That's quite wrong - MT requires a box which must be AT LEAST double the size of the electronic density.
Providing only 3-5 or 6 Angstrom more for each dimension will lead to wrong results. Or would work for molecule of maximum 10 Angstrom wide (in the direction in which you add that buffer).
Other methods, not implemented in CP2K, like Hockney require that the density has to be zero at the
edges of the box. This is not the case for MT: it is even more demanding.
So if the molecule is 15 Angstrom in one dimension (I'm talking only about atomic distances) the minimum box will be 30 + some buffer (3-5 Angstrom) to keep into account the electronic density in that dimension.
Also using an adaptive box can be quite dangerous in terms of convergence of results: if the molecule is flat one has to consider only the thickness of the electronic density and without knowing it a priori any number (3 or 5) can be absolutely arbitrary.
Best,
Teo
> To unsubscribe from this group, send email to cp2k+uns...@googlegroups.com.
Matthias Krack
you are right, I was mistaken, mixing Hockney and MT.
sorry,
Matthias
>
> read more »
Teodoro Laino
the script you've been using should be the one we distributed at the CP2K tutorial.
It's strange that Matthias does not know anything about that since he was tutor at that tutorial.
Anyway: the script does exactly what I said in the previous mail.
It takes the linear dimensions of the molecule (atomic distances) and sums up a buffer 3-5 Angstrom.
This will lead to converged results for decouplers like MT.
Again in CP2K we do not have cheaper decoupler (you may try the others one available like :
MULTIPOLE
WAVELET
)
Refer to the papers cited for each of these methods for constraints about the box size: each of them is slightly different.
For large molecules WAVELET may be definitely much better. In fact for WAVELET is also enough that the density is zero at the border of the box.
But there are other constraints : for non-periodic calculations we have implemented only cubic boxes. So.. you see.. depends a lot..
for MT you can use an orthorhombic (but larger in size).. for WAVELET a cubic although a bit smaller.
Up to you the final decision.
Teo
Jörg Saßmannshausen
many thanks for your replies.
Yes, I remember from the tutorial that you need the 'right size of the box'
and, as mentioned by Theo, should be at least double the size of the electron
density. Fortunately, I should have access to a larger cluster over Christmas.
Now, correct me if I am wrong, the memory of the whole calculation should be
(fairly) the same, say 64 GB, regardless whether I run it hypothetically on
one or 64 cores, right? I am aware of some overhead when I run things in
parallel. The reason I am asking is, the IB cluster here has 4 cores per node
with 8 GB or RAM, so I would need some rule of thumb to estimate the numbers
of cores I would need.
The reason behind all of that is that I need to do some testing to get the
right hardware configuration for the new cluster we want to purchase.
Obviously, I want enough memory (I am currently applying 3 GB per core for the
Intel Nehalem) but I don't want to overkill it with memory. Hence all my
testing.
All the best from a cold London
Jörg
> >> ON. html to get an idea of what will be used.
Matt W
It is definitely worth looking closely at the RS grids. One (more)
technicality is that, by default the 2 highest RS grids are
distributed across the processors, thus total memory should conserved,
as it were and reduce per core with #procs. Grids below this, I
think, are fully replicated across all processors, which might not be
a good thing for your big box - it will give a large memory term that
will not reduce with #procs. This can be altered in the MGRID
section, but hard to say more without actually seeing data.
Best from rather grim London,
Matt
On Dec 17, 9:56 am, Jörg Saßmannshausen <j.sassmannshau...@ucl.ac.uk>
>
> read more »
