Dear Alessandro,
thank you for this explanation!
Kind regards,
Georg

Hi Alessandro,
thank you very much for these detailed explanations!
Crystal is an amazing program, and I am just starting to explore its powerful thermoelastic capabilities.
Kind regards,
Georg

Hi Erba,
Honestly, I have ANDERSON acceleration on because I inherited a file with these keywords written in when I was taught how to use CRYSTAL for frequency calculations - and worked pretty well so far with our other systems so I thought to leave it alone (though it was more for pharmacuetical systems). Thanks for the suggestion, I'll try it out if with the original TOLINTEG 10 10 10 10 20 and see what happens.

Thanks,
Peter

Alessandro, the setback is temporary. What I really appreciate is that you gave CRYSTAL more life via user support than it ever had. I am very grateful for it

Hi,

The first data series corresponds to the first set and the second to the second.
If you need more insight on the specifics of your calculation, please share your files.

Hope this helps,

Hi,

In CRYSTAL, the most general syntax to specify exchange-correlation (xc) functionals is, within the DFT input block through the EXCHANGE and CORRELAT keywords as:

DFT EXCHANGE label of exchange functional (e.g. VBH, BECKE, PBE, ...) CORRELAT label of correlation functional (e.g. VWN, LYP, PBE, ...) ENDDFT

For certain standard combinations of exchange and correlation functionals, we have implemented single keywords. For instance:

DFT BLYP ENDDFT

is equivalent to

DFT EXCHANGE BECKE CORRELAT LYP ENDDFT

Another example (here the "XC" letters are appended at the end of the name to reflect that the functional is used for both the exchange and correlation part):

DFT PBESOLXC ENDDFT

is equivalent to

DFT EXCHANGE PBESOL CORRELAT PBESOL ENDDFT

For PBE, there are three equivalent ways to define it in CRYSTAL:

DFT PBEXC ENDDFT

or

DFT PBE ENDDFT

or

DFT EXCHANGE PBE CORRELAT PBE ENDDFT

There isn't a general rule, I'm afraid. This syntax aspects are discussed at page 134 of CRYSTAL23 User's Manual.

Hope this clarifies things a little,

Hi,

Direct piezoelectric constants of a 3D lattice in CRYSTAL are defined and computed as:
$$
e_{ci}^{3D} = \left( \frac{\partial P_c}{\partial \eta_i}\right) = \frac{1}{V}\left( \frac{\partial^2 E}{\partial E_c\partial \eta_i}\right)
$$
that is as first derivatives of Cartesian components of the polarization (c=x,y,z) with respect to strain components, or, equivalently as second derivatives of the energy density (V is the volume of the 3D lattice cell) with respect to Cartesian components of an electric field \(E_c\) and strain components, where the strain \( \eta \) is dimensionless and thus the direct piezoelectric constants have units of \( \textup{charge/length}^2 \).

For 1D and 2D periodic lattices, as the volume (V) is not uniquely defined (or not defined at all in some cases), one may divide by the length \(l \) and area \( A\) of the lattice cell instead:
$$
e_{ci}^{1D} = \frac{1}{l}\left( \frac{\partial^2 E}{\partial E_c\partial \eta_i}\right) \quad \textup{and} \quad e_{ci}^{2D} = \frac{1}{A}\left( \frac{\partial^2 E}{\partial E_c\partial \eta_i}\right)
$$
that would thus be expressed in units of \( \textup{charge} \) or \( \textup{charge/length} \) for 1D and 2D lattices, respectively.

However, in CRYSTAL for 1D and 2D lattices we do not divide by \(l \) or \( A\) , and just define and compute the piezoelectric constants as:
$$
e_{ci}^\textup{1D and 2D} = \left( \frac{\partial^2 E}{\partial E_c\partial \eta_i}\right)
$$
with units of \( \textup{charge}\cdot\textup{length} \).

Yes, these constants are physically meaningful for 1D and 2D systems. For a 2D monolayer system, for instance, depending on what you need to compare with, you can do one of two things:

keep them as they are printed in the CRYSTAL output (units of \( \textup{charge}\cdot\textup{length} \))

divide the values you get in the CRYSTAL output by the area of the 2D cell (and thus express them in units of \( \textup{charge/length} \))

I would not divide by a volume because I would not know the physical meaning of the volume of a 2D monolayer system.

Hope this helps,

Hi,

The MP2 option is no longer supported in recent versions of the CRYSTAl program. If you are interested in a periodic MP2 calculation, my suggestion is to contact Lorenzo Maschio ([email protected]) and Denis Usvyat ([email protected]) directly, who may provide guidance in the use of the CRYSCOR program.

I tried rerunning it with fewer nodes - thought it is some parallel issue. A problem again

ANGULAR INTEGRATION - INTERVALS (ACCURACY LEVEL [N. POINTS] UPPER LIMIT):
1( 4[ 86] 0.2) 2( 8[ 194] 0.5) 3( 12[ 350] 0.9) 4( 16[ 974] 3.5)
5( 12[ 350]9999.0)
CYCLE 0 ALPHA 227.814788 EPSILON 1.894274 DELTA 2.2781E+02
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 1902.88 TCPU 1885.42
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 1904.27 TCPU 1886.79
DIIS TEST: 0.61205E+01 AT CPHF CYCLE 1 - MIX 60 %
CYCLE 1 ALPHA 257.133404 EPSILON 2.009363 DELTA 2.9319E+01
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2002.77 TCPU 1984.78
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2004.16 TCPU 1986.16
DIIS TEST: 0.71887E+01 AT CPHF CYCLE 2 - DIIS ACTIVE - HISTORY: 2 CYCLES
CYCLE 2 ALPHA 268.265588 EPSILON 2.053062 DELTA 1.1132E+01
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2102.04 TCPU 2083.54
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2103.42 TCPU 2084.92
DIIS TEST: 0.36370E+00 AT CPHF CYCLE 3 - DIIS ACTIVE - HISTORY: 3 CYCLES
CYCLE 3 ALPHA 276.769385 EPSILON 2.086443 DELTA 8.5038E+00
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2202.03 TCPU 2183.04
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2203.42 TCPU 2184.42
DIIS TEST: 0.54051E-01 AT CPHF CYCLE 4 - DIIS ACTIVE - HISTORY: 4 CYCLES
CYCLE 4 ALPHA 278.095061 EPSILON 2.091647 DELTA 1.3257E+00
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2302.12 TCPU 2282.64
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2303.51 TCPU 2284.02
DIIS TEST: 0.85023E-02 AT CPHF CYCLE 5 - DIIS ACTIVE - HISTORY: 5 CYCLES
CYCLE 5 ALPHA 278.435921 EPSILON 2.092985 DELTA 3.4086E-01
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2402.16 TCPU 2382.20
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2403.54 TCPU 2383.57
DIIS TEST: 0.38480E-03 AT CPHF CYCLE 6 - DIIS ACTIVE - HISTORY: 6 CYCLES
CYCLE 6 ALPHA 278.461661 EPSILON 2.093086 DELTA 2.5739E-02
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2502.06 TCPU 2481.62
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2503.45 TCPU 2482.99
DIIS TEST: 0.44991E-03 AT CPHF CYCLE 7 - DIIS ACTIVE - HISTORY: 7 CYCLES
CYCLE 7 ALPHA 278.460154 EPSILON 2.093080 DELTA -1.5071E-03
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2601.70 TCPU 2580.74
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2603.08 TCPU 2582.11
DIIS TEST: 0.36243E-03 AT CPHF CYCLE 8 - DIIS ACTIVE - HISTORY: 8 CYCLES
CYCLE 8 ALPHA 278.473843 EPSILON 2.093134 DELTA 1.3689E-02
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2701.77 TCPU 2680.26
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2703.15 TCPU 2681.62
DIIS TEST: 0.85073E-04 AT CPHF CYCLE 9 - DIIS ACTIVE - HISTORY: 9 CYCLES
CYCLE 9 ALPHA 278.474328 EPSILON 2.093136 DELTA 4.8487E-04
forrtl: severe (256): unformatted I/O to unit open for formatted transfers, unit 85, file /dev/null
Image PC Routine Line Source
Pcrystal 0000000007374206 Unknown Unknown Unknown
Pcrystal 0000000001BA179E Unknown Unknown Unknown
Pcrystal 0000000000A8038B Unknown Unknown Unknown
Pcrystal 0000000000A63D97 Unknown Unknown Unknown
Pcrystal 0000000000D4DAD1 Unknown Unknown Unknown
Pcrystal 000000000074B942 Unknown Unknown Unknown
Pcrystal 000000000040591E Unknown Unknown Unknown
Pcrystal 00000000004053FD Unknown Unknown Unknown
libc.so.6 000014B7C14295D0 Unknown Unknown Unknown
libc.so.6 000014B7C1429680 __libc_start_main Unknown Unknown
Pcrystal 0000000000405315 Unknown Unknown Unknown

Hello Dr. Erba,

Thank you very much for your thorough explanation! This clarifies all my questions. I will make sure to re-calculate using a larger supercell.

Best,
Danny

Hi!

In a COOP calculation you aim to have a description on the interaction between pairs of orbitals or atoms. In order to do that, you need to indicate each pair you are interested to analyze. In your input, in the first line after the COOP keyword, the initial number 1 indicates that you are interested in one pair of orbitals/atoms. You still need to indicate a pair of orbitals or atoms to be considered, writing them in separated lines. Consider this example, taken from the Tutorials webpage:

NEWK 6 6 1 0 COOP 1 200 7 14 1 12 0 -1 1 -1 2 END

Here, the two lines before the final END keyword indicate which atoms will be considered (atoms, given that the lines start with a negative value, as stated in the manual page 322). COOP will be evaluated considering the first and second atoms of the systems (with indices 1 and 2). From your previous calculations you can recover the indices of the atoms/orbitals you are interested.

Let me know if this information has been useful 🙂

Thank you very much. Now it works without any issue. I was not aware of the necessity to change this parameter as I did not encounter any issue with other calculations I have done due to it being 0 0 0.

OK, here we go. It just is stuck, always the same position in the output

(ceres20-compute-46:0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95)
(ceres24-compute-18:96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191)

export TMPDIR=/local/bgfs/jonas.baltrusaitis/15383115
export TMOUT=5400
export SINGULARITY_TMPDIR=/local/bgfs/jonas.baltrusaitis/15383115

MAX NUMBER OF SCF CYCLES 200 CONVERGENCE ON DELTAP 10**-20
WEIGHT OF F(I) IN F(I+1) 30% CONVERGENCE ON ENERGY 10**-10
SHRINK. FACT.(MONKH.) 6 6 6 NUMBER OF K POINTS IN THE IBZ 64
SHRINKING FACTOR(GILAT NET) 6 NUMBER OF K POINTS(GILAT NET) 64

*** K POINTS COORDINATES (OBLIQUE COORDINATES IN UNITS OF IS = 6)
1-R( 0 0 0) 2-C( 1 0 0) 3-C( 2 0 0) 4-R( 3 0 0)
5-C( 0 1 0) 6-C( 1 1 0) 7-C( 2 1 0) 8-C( 3 1 0)
9-C( 0 2 0) 10-C( 1 2 0) 11-C( 2 2 0) 12-C( 3 2 0)
13-R( 0 3 0) 14-C( 1 3 0) 15-C( 2 3 0) 16-R( 3 3 0)
17-C( 0 0 1) 18-C( 1 0 1) 19-C( 2 0 1) 20-C( 3 0 1)
21-C( 0 1 1) 22-C( 1 1 1) 23-C( 2 1 1) 24-C( 3 1 1)
25-C( 0 2 1) 26-C( 1 2 1) 27-C( 2 2 1) 28-C( 3 2 1)
29-C( 0 3 1) 30-C( 1 3 1) 31-C( 2 3 1) 32-C( 3 3 1)
33-C( 0 0 2) 34-C( 1 0 2) 35-C( 2 0 2) 36-C( 3 0 2)
37-C( 0 1 2) 38-C( 1 1 2) 39-C( 2 1 2) 40-C( 3 1 2)
41-C( 0 2 2) 42-C( 1 2 2) 43-C( 2 2 2) 44-C( 3 2 2)
45-C( 0 3 2) 46-C( 1 3 2) 47-C( 2 3 2) 48-C( 3 3 2)
49-R( 0 0 3) 50-C( 1 0 3) 51-C( 2 0 3) 52-R( 3 0 3)
53-C( 0 1 3) 54-C( 1 1 3) 55-C( 2 1 3) 56-C( 3 1 3)
57-C( 0 2 3) 58-C( 1 2 3) 59-C( 2 2 3) 60-C( 3 2 3)
61-R( 0 3 3) 62-C( 1 3 3) 63-C( 2 3 3) 64-R( 3 3 3)

DIRECT LATTICE VECTORS COMPON. (A.U.) RECIP. LATTICE VECTORS COMPON. (A.U.)
X Y Z X Y Z
13.1430453 0.0000000 0.0000000 0.4780616 0.0000000 0.0000000
0.0000000 11.6066979 0.0000000 0.0000000 0.5413413 0.0000000
0.0000000 0.0000000 21.1989478 0.0000000 0.0000000 0.2963914

DISK SPACE FOR EIGENVECTORS (FTN 10) 53868000 REALS

SYMMETRY ADAPTION OF THE BLOCH FUNCTIONS ENABLED
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT gordsh1 TELAPSE 186.18 TCPU 45.44

Hi Othmen!

I think that you have a problem in the geometry sections of your input files. When dealing with monoclinic structures, only one angle (beta) has to be specified as alfa and gamma are 90 by default, see page 22 of Crystal User's Manual. Furthermore, in input1 the atomic coordinates are missing, while in input2 94 atoms coordinates are written but you specified that 95 atoms are present.

Let me know if you have other problems 🙂

For the purpose of finding the minimum energy structure to then do Raman calculations, it is.

EOS gives you much more than that of course: the p(V) or, equivalently, V(p) relation (i.e. structure as a function of pressure), the bulk modulus K(p), and allows to compute the enthalpy H(p).

Hi aimipa!

Yes there is an option to set an anisotropic grid from NEWK, you can find it at page 346 of the CRYSTAL User's Manual, see also the screenshot I attach below. You have to set IS=0 and then you can choose three different shrinking factors along B1, B2 and B3.

For example:

NEWK 0 12 12 8 4

sets a shrinking factor of 12, 8 and 4 along the three reciprocal lattice vectors.

I have never tested this anisotropic option in a BOLTZTRA calculation, would you mind let me know if it works? 🙂

newk.png

Good! The clean way is:

CRYSTAL
1 0 0
P 21/A

I guess the extra "1"s you put are safe as they correspond to the identity operator.

Yes sure. No problem.
I'm copying it here below.
Thank you for the help

Daria

Furfural CRYSTAL 0 0 0 19 7.83542762 10.92181175 14.94182734 90.000000 90.000000 90.000000 33 8 1.907561246186E-01 3.035506248522E-01 4.408415126674E-01 8 -4.953260939445E-01 3.221701934644E-01 3.435788389827E-01 6 2.036695655476E-01 3.000622663369E-01 3.482645168957E-01 6 4.363374159111E-02 2.908170991680E-01 3.122160878064E-01 1 1.581732742072E-02 2.880586476013E-01 2.419150131741E-01 6 -7.402527715589E-02 2.893444053747E-01 3.842913887226E-01 1 -2.109486401803E-01 2.822752418693E-01 3.798189138779E-01 6 2.099722211847E-02 2.968688203007E-01 4.606664256414E-01 1 -9.310483322231E-03 2.937832972654E-01 -4.691825366176E-01 6 3.647591498990E-01 3.132132214593E-01 3.052219591009E-01 1 3.550467304924E-01 3.168499490771E-01 2.320173417360E-01 8 2.803460822693E-01 -3.091160241316E-01 4.442088720927E-01 8 -4.110526913885E-01 -3.641439806750E-01 3.516630126321E-01 6 2.872037985115E-01 -3.695128622796E-01 3.626185432797E-01 6 1.245767077225E-01 -3.991101749711E-01 3.362778399038E-01 1 9.163169834882E-02 -4.467733390550E-01 2.755816936910E-01 6 1.226510841197E-02 -3.550438914407E-01 4.034527830796E-01 1 -1.250093786328E-01 -3.609739865930E-01 4.033488316714E-01 6 1.122081912874E-01 -3.009269785992E-01 4.671425614917E-01 1 8.596288473646E-02 -2.539158941058E-01 -4.709004162095E-01 6 4.474098992917E-01 -3.907868742186E-01 3.203804132925E-01 1 4.349587971852E-01 -4.336723578587E-01 2.540887571847E-01 8 1.068738713332E-01 -3.851604875541E-02 3.194441316968E-01 8 -2.354778627755E-01 -1.170305051641E-02 3.709929081706E-01 6 5.921416878240E-02 -5.109842318475E-03 4.052174692072E-01 6 2.027905991547E-01 1.001107546720E-02 4.563966208009E-01 1 2.038099671676E-01 3.508046669402E-02 -4.738476409789E-01 6 3.453540637796E-01 -1.523126326522E-02 4.003413693017E-01 1 4.790458182320E-01 -1.335338648893E-02 4.173410231626E-01 6 2.808235437569E-01 -4.461452496729E-02 3.184022003331E-01 1 3.367732131084E-01 -7.306049141011E-02 2.558072164743E-01 6 -1.187612609562E-01 6.696688015001E-03 4.250681762815E-01 1 -1.454796120848E-01 3.516201574611E-02 4.940534466886E-01 FREQCALC NUMDERIV 2 INTENS INTRAMAN INTCPHF END END BASISSET POB-TZVP-REV2 DFT B3LYP-D3 XLGRID ENDdft SHRINK 4 4 TOLDEE 10 SCFDIR LEVSHIFT 6 1 FMIXING 30 TOLINTEG 7 7 7 7 14 EXCHSIZE 13220300 BIPOSIZE 13220300 ENDscf

Hi,

Hmmm...even though the SI file of the paper you have shared does point to using HSE06, it remains not entirely clear whether its parameters were kept default or altered (for all systems together or individually for instance).

There are also differences between what different codes consider as "default" settings. For example, within the code used in the paper (VASP, nice code, no doubt), the defaults for HSE06 read (taken from https://www.vasp.at/wiki/index.php/List_of_hybrid_functionals) :
$$ \omega= 0.2\ \mathring{A} , \quad c = 0.25, \quad \text{correlation}=\text{PBE}, $$ with the first number reading the range separation parameter (omega) and the second the fraction of exact exchange used (c).

Within CRYSTAL (also nice code, no doubt), these read (taken from the manual, page 138):
$$ \omega= 0.11\ a_0^{-1}, \quad c = 0.25, \quad \text{correlation}=\text{PBE}, $$ adopting the same labels.

Not sure about the exact definition of units (perhaps a developer can comment if this is indeed Bohr radius as assumed?), but you can already see the subtle differences having to be taken into account when comparing between codes.

A few other thoughts worth considering:

In the paper, the structure was optimized with PBEsol and on top of that geometry HSE06 was applied as a single-point calculation. Not sure about the exact composition of those ZIFs, but the structural differences could play a significant role as well (planewave codes are very costly when optimizing a structure with hybrid functionals). Here is also a good read on this topic: doi.org/10.1088/2516-1075/aafc4b

One final small comment. Within the PAW formalism implemented in VASP, scalar relativistic effects are included in the pseudopotentials by default. No problem, cool feature, but should be taken into account when comparing results, especially for heavier elements (longer discussion found here https://blog.vasp.at/forum/viewtopic.php?t=902)

Hope this helps!

Cheers,
Aleks

Thanks Jacques! These are very informative.
The information for specific basis sets is helpful, I will try them.