Yes, indeed, the d-type shell is missing. This is so unfortunate....

Worst case scenario is the only possible scenario it seems...

Although it is difficult to predict what effect this missing d orbitals will have on the computed frequencies, it may be significant.

Ok, so, you'll need to insert the basis set explicitly in the input for all elements in your system. That is, in the standard way, without the BASISSET keyword. For S, the complete pob-tzvp-rev2 basis is:

16 10
0 0 7 2.0 1.0
60700.928104 0.00054695944225
9102.6106854 0.00422972245570
2071.4166009 0.02174782415900
586.02476821 0.08510005358900
190.55395021 0.24799128459000
67.630384260 0.46703640406000
25.127306905 0.36434587550000
0 0 3 2.0 1.0
112.57463010 0.02167004024000
34.795554217 0.09360230176000
6.5115556215 -0.26068001422000
0 0 2 2.0 1.0
3.2399032261 1.28420894350000
1.5477160881 0.66036416584000
0 0 1 0.0 1.0
0.4487335200 1.00000000000000
0 0 1 0.0 1.0
0.1553457200 1.00000000000000
0 2 5 6.0 1.0
564.36716027 0.00247967963170
133.42624379 0.01967793025000
42.468271189 0.08998000825800
15.616527580 0.25705880575000
6.1093988469 0.43515167292000
0 2 1 4.0 1.0
2.0359436000 1.00000000000000
0 2 1 0.0 1.0
0.4337928300 1.00000000000000
0 2 1 0.0 1.0
0.1305009100 1.00000000000000
0 3 1 0.0 1.0
0.4107010100 1.0000000000000

You can find the other basis sets and copy-paste them here

Although this is not my personal fault, I feel very sorry for this ugly bug.

Hi,

Before sharing my general thoughts about POB versus "old school" basis sets, let me ask if you are using POB-TZVP or POB-TZVP-REV2 for your Sulfur-containing systems.

It was brought to my attention just a couple of days ago that the POB-TZVP-REV2 basis for Sulfur (S) is missing a whole d-type shell in the internal libraries in CRYSTAL and is thus not to be trusted, unless explicitly inserted in the input file as copy-pasted from the website!!!

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,

By the looks of it, it may be due to some mismatch in the restart units. I am happy to run some tests. Could you share your input files for this case?

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,

I have run some calculations on your case following my step-by step recipe (as described at https://forum.crystalsolutions.eu/post/339, which I have now edited to include the intensity step), and the restart now works just fine (i.e. without the annoying "possibly conducting state", and with a nice convergence of the further SCFs upon restart).

This is what I did:

• I took the optimized geometry from your original output file and I created a new input file for a single-point calculation. I ran it and obtained the wavefunction external file (i.e. fort.9 unit);

• I prepared an input file to run the harmonic frequency calculation, I ran it and I killed it in the middle of the construction of the Hessian. From this incomplete frequency calculation, I obtained the FREQINFO.DAT file and the external unit with the density matrix, i.e. fort.13 unit;

• I prepared an input file to restart the frequency calculation and provided the FREQINFO.DAT and fort.13 files from the previous step and the fort.9 (renamed as fort.20) from the first step, and ran it. The calculation of the Hessian restarted correctly, with the new SCFs converging nicely (see the attached SCFOUT.LOG file). I then stopped the calculation not to use too much compute power on my cluster.

If you follow this step-by-step process you'll be able to safely restart your frequency calculations.

Hope this helps,

job314 I am running some tests now. I reproduced your "possibly conducting state" problem. Let me try a couple of things.

Yes, please, send me all the files.

Hi,

I am quite familiar with that error message myself: it pops up when the fort.9 unit provided upon the restart is the wrong one (i.e. does not match with the current calculation).

Let me make a general comment on how to best approach these frequency+intensity calculations with CRYSTAL based on my experience.

I noticed that you tend to use the PREOPTGEOM option within FREQCALC, and to compute infrared and Raman intensities, all in one shot. Now, this is a very nice feature of CRYSTAL that allows you to run a single job where everything is fully automated: the structure gets optimized, the numerical Hessian computed and diagonalized, and Born and Raman tensors computed. Not many programs can do that, to the best of my knowledge. At the same time, while this is very convenient if you can indeed run everything in a single job, it complicates things if you then need to do a restart from a previous incomplete calculation.

If you envisage that this could be the case (maybe because of a wall clock limit on the cluster) then it is preferable (to put it mildly) to do things step by step:

• You first run a geometry optimization with something like:

[initial geometry]

OPTGEOM
END

• You prepare a new input file where you insert the optimized geometry from the previous run as a starting point, and perform a frequency calculation, with something like:

[optimized geometry]

FREQCALC
END

• If this calculation stops, you can now restart it easily with:

[optimized geometry]

FREQCALC
RESTART
END

by providing the required restart files (FREQINFO.DAT and fort.13 from the first frequency calculation, and fort.9 from the geometry optimization, to be renamed fort.20 in the new scratch folder).

• Once the restart is done and the frequency calculation is complete, you can compute the intensities, with something like:

[optimized geometry]

FREQCALC
RESTART
INTENS
INTRAMAN
INTCPHF
ENDCPHF
END

Personally, I always do things separately, running first the geometry optimization, then the Hessian (i.e. harmonic frequency) calculation, then the intensities through CPHF/KS in three separate jobs.

Hope this clarifies things a little,

Whatever works, I guess

The point that you make about the geometry may also be relevant. Given that your original calculation was also performing a geometry optimization, when you do the restart calculation you must ensure that the geometry specified in the new input is the optimized one from the previous run.

Hi,

Let's try to fix one of these restart issues at a time.

As you know (I put it here for other readers as well), to restart a FREQCALC calculation we need the RESTART keyword, as below:

FREQCALC
RESTART
END

All keywords related to the calculation of the intensities from the initial job should be kept in the new input as well.

In order to make this restart work, you need a few files from the previous job to be placed in the scratch folder of the new job: FREQINFO.DAT, fort.13 and fort.9 (to be renamed fort.20 in the new folder).

As such, I am unable to try to restart your job because I need these three files from the previous run. If you have them you can try to restart and/or send me the files.

Let's see how this goes.

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.

Hi,

Let me try to rationalize what you observe as follows: In CRYSTAL, interatomic force constants for phonon dispersion are computed with a direct space approach within a supercell. Because of the reciprocity relation between direct and reciprocal lattices, the number of k-points for which phonons can be computed is determined by the expansion factor of the supercell relative to the primitive cell. For instance, with your 2x2 supercell expansion for a 2D system, you get to compute phonons at just 4 k-points. In order to get the phonons at more k-points you can do one of two things (or both):

1. Increase the supercell expansion (i.e. increase the size of the supercell). Going to 3x3 gives you phonons at 9 k-points, 4x4 gives you phonons at 16 k-points, and so on. This is done through the SCELPHONO keyword:

SCELPHONO
6 0
0 6

2. Interpolate the interatomic force constants to evaluate the dynamical matrices at additional k-points with respect to those determined by the expansion matrix of the supercell. This can be done through the INTERPHESS option as follows:

INTERPHESS
8 8
0

NOTE 1: This interpolation is safe only if the interatomic force constants vanish within the supercell. This usually requires the supercell to have a radius of at least 10-15 Angstroms (i.e. the supercell has to be large).

NOTE 2: This interpolation is used when phonon bands are computed with BANDS. Thus, the starting supercell has to be large to compute phonon bands safely.

In your case, what you observe is the following: when you do the BANDS calculation you are doing the interpolation. When you do the PDOS calculation you are not. To add the interpolation in your PDOS calculation (and get a phonon density-of-states that matches that of the phonon band plot) you need to add the INTERPHESS option.

But, and this is an important but, in this way you would get matching results but they are likely to be both wrong because the supercell you are starting from is too small. To get more reliable results, you need to significantly increase the supercell expansion matrix.

Hope this helps,

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).

Absolutely yes, the optimized structure you get from OPTGEOM (i.e. a full optimization of both atomic coordinates and cell shape/volume, in CRYSTAL) is very very very close the minimum from fitted EOS, and if they differ I would blame EOS over OPTGEOM (meaning that maybe more points in the EOS scan would be required). Additionally, OPTGEOM is much much faster than EOS as only 1 versus N optimization processes are actually performed.

So, to get the fully relaxed structure, I strongly recommend to use OPTGEOM in CRYSTAL. This is what I usually do before running frequency and Raman calculations.

Maybe I understand where the confusion may come from. Till CRYSTAL09, the OPTGEOM option corresponded to a geometry optimization of just the atomic coordinates within a fixed cell. Since CRYSTAL14 (and thus in versions 14, 17 and 23) the OPTGEOM option corresponds to a full geometry optimization (i.e. OPTGEOM is now equivalent to FULLOPTG).

So, OPTGEOM provides the minimum energy structure.

I would need to better understand what you mean by systematic scanning of lattice parameters to comment on that point.

Let me just reiterate on the numerical nature of both the OPTGEOM and EOS options: the minimum structure from both approaches should be the same. If they differ, it is just due to numerical reasons and, personally, I would tend to trust OPTGEOM better.

Hi,

I think that my point is that if you want to find the minimum of the PES you really do not need to run an EOS calculation. You can simply run an OPTGEOM calculation. EOS is useful if you need to take the pressure into explicit account and/or compute the bulk modulus.

If you have already run an EOS calculation with the PREOPTGEOM option, the structure corresponding to the minimum is the one obtained at the end of the pre-optimization, thus before all the various 0.94, 0.96, etc. compressions/expansions are explored. So in this case, all you need to do is look in your output file for the first occurrence of the "OPT END" string. You will find the optimized structure printed in its vicinity.

Hope this helps,

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.