Hi,

Jefferson Maul has managed to generate this .cif file with 4 symmetry operators. As you suggest, there are probably more but this is what he could extract so far.

Beautiful system by the way: looks like a Christmas tree bauble!

Hi,

job314 said in SCANMODE io error Read_int_1d:

So it turns out having fort.13 and fort.20 is not optional for restart.

Yes, the lack of these two files was the origin of the I/O error.

job314 said in SCANMODE io error Read_int_1d:

PS. I would also like to ask developers to allow uploading compressed files to this forum

Now also .zip, .tar, .tgz, .tar.gz files can be uploaded.

Cheers,

Hi,

Before running the actual single-point or geometry optimization calculations on the configurations with RUNCONFS, a list of configurations needs to be generated by use of the CONFRAND option. The list of generated configurations is saved into a file CONFIGURATIONS.DAT that is then read by the next RUNCONFS calculation.

Let us go through this step-by-step. I take your system as an example.

• First, you would setup an input for the CONFRAND calculation. For instance:

Title
CRYSTAL
0 0 0
194
3.065 17.656
4
22 0 0 0.5
14 0 0 0.75
22 0.666666 0.333333 0.364919
6 0.333333 0.666666 0.427507
SCELCONF
1 0 0
0 1 0
0 0 1
CONFRAND
1
5
2
END

Here I am selecting just one crystallographic site for substitution, specified by atom number 5, which is a Ti atom. Given the symmetry of this system, that atom has a multiplicity of 4 (i.e. there are other 3 Ti atoms symmetry-related to it). This can be inspected from here (in bold the selected atoms, in italic its symmetry-equivalents):

 N. ATOM EQUIV AT. N.          X                  Y                  Z

   1   1   1   22 TI    0.00000000000E+00  0.00000000000E+00 -5.00000000000E-01
   2   1   2   22 TI    0.00000000000E+00  0.00000000000E+00  0.00000000000E+00

   3   2   1   14 SI    0.00000000000E+00  0.00000000000E+00 -2.50000000000E-01
   4   2   2   14 SI    0.00000000000E+00  0.00000000000E+00  2.50000000000E-01

 **5   3   1   22 TI   -3.33334000000E-01  3.33333000000E-01  3.64919000000E-01**
  *6   3   2   22 TI    3.33334000000E-01 -3.33333000000E-01 -1.35081000000E-01
   7   3   3   22 TI    3.33333000000E-01 -3.33333000000E-01 -3.64919000000E-01
   8   3   4   22 TI   -3.33333000000E-01  3.33334000000E-01  1.35081000000E-01*

   9   4   1    6 C     3.33333000000E-01 -3.33334000000E-01  4.27507000000E-01
  10   4   2    6 C    -3.33333000000E-01  3.33334000000E-01 -7.24930000000E-02
  11   4   3    6 C    -3.33333000000E-01  3.33334000000E-01 -4.27507000000E-01
  12   4   4    6 C     3.33334000000E-01 -3.33333000000E-01  7.24930000000E-02

The last input parameter of CONFRAND, which I set to 2, determines how many of these 4 Ti atoms will be substituted.

• By running it you get the following output:

*******************************************************************************
  SUBSTITUTIONS AT SITES (LABELS) : 
    5    7    6    8
  **********************************            COMPOSITION :          2 /   4
  **********************************          NUMBER OF SIC :                3
 *******************************************************************************
 --->              1      SIC  FOUND AT TRY              1   -   CONFIGURATION 
   0   0   1   1
   MULTIPLICITY      2   -   RANK           1   -   CANONICAL RANK           1
 --->              2      SIC  FOUND AT TRY              2   -   CONFIGURATION 
   1   0   1   0
   MULTIPLICITY      2   -   RANK           5   -   CANONICAL RANK           2
 --->              3      SIC  FOUND AT TRY              4   -   CONFIGURATION 
   1   0   0   1
   MULTIPLICITY      2   -   RANK           3   -   CANONICAL RANK           4
 *******************************************************************************
                   3      SIC  FOUND  AFTER              4       TRIES
 *******************************************************************************

that is 3 symmetry-independent configurations (SICs) are found and stored in the external file CONFIGURATIONS.DAT.

• At this point you are ready to run a RUNCONFS calculation (note that the CONFIGURATIONS.DAT file generated at the previous step needs to be placed inside the scratch folder of the new job). For instance with:

Title
CRYSTAL
0 0 0
194
3.065 17.656
4
22 0 0 0.5
14 0 0 0.75
22 0.666666 0.333333 0.364919
6 0.333333 0.666666 0.427507
SCELCONF
1 0 0
0 1 0
0 0 1
RUNCONFS
ATOMSUBS
22 273
END

In this case I ask to substitute Ti with Ta.

Hope this helps,

Hi,

Can you check the COHP.dat file you uploaded? I tried to open it with CRYSPLOT but everything seems to be zero...?

Hi,

jquertin said in SCF fails spinlock with POB-DZVP-REV2:

In the case of SPINLOCK, the manual clearly explain that the NSPIN value is the difference in number of alpha and beta electrons. For SPINLOC2, the text only refers to the spin while the table gives the same definition for SPIN as for NSPIN (in SPINLOCK).

The argument SPIN of SPINLOC2 still represents a number of electrons, as in SPINLOCK.

jquertin said in SCF fails spinlock with POB-DZVP-REV2:

Furthermore, in the calculation, if using SPINLOC2 with 6 or 6.0 as the spin (as defined in the table), crystal defaults to SPINLOCK.

That's right. SPINLOC2 requires a non integer argument. For integer arguments it reduces to SPINLOCK.

jquertin said in SCF fails spinlock with POB-DZVP-REV2:

In short, if I define SPINLOC2 SPIN as 3 (1/2 * 6) or 3.0, crystal defaults to SPINLOCK with NSPIN 3 which is actually half of what I want.

In both SPINLOCK and SPINLOC2, the argument is meant as a number of electrons. Thus, if you have 6 extra up electrons with respect to down electrons, the input value should be 6, not 3. For integer values, SPINLOC2 is of no use.

Hope this clarifies things a little,

Hi,

To specify your plane of interest more freely, I recommend using the COORDINA option instead of the ATOMS option. With COORDINA you can input the Cartesian coordinates of 3 points A, B and C defining the plane. For instance, if you want to select the XY plane you can simply use:

ECHG
0
50
COORDINA
1.0 0.0 0.0
0.0 0.0 0.0
0.0 1.0 0.0
RECTANGU
MARGINS
3 3 3 3
END
END

Hope this helps,

Hi,

Just a couple of comments:

1. Are you starting from a fully optimized structure? I ask because the following warning is printed in your output file:

WARNING!! FORCE   9 AT CENTRAL POINT IS GREATER THAN 10^-4. OPTIMIZE THE STRUCTURE AND RE-RUN

2. Anisotropic shrinking factors are not fully supported by ELAPIEZO calculations, which may be at the origin of the problem you are experiencing. I would advise to try running the calculation with an isotropic shrinking factor. Something like:

SHRINK
3 3

Let us know if this helps,

Hi,

1. There may not be an option to "plot" orbital-resolved bands directly, but there is an option to compute and print principal atomic orbitals contributions to selected eigenvectors. See the ANBD keyword of the PROPERTIES module (page 307 of the CRYSTAL23 User's Manual).

2. Yes, if your system is 2D, you can use the usual BAND keyword of the PROPERTIES module to define a 2D path to compute and plot the band structure (page 309 of the CRYSTAL23 User's Manual).

If you need help on a specific system, just let me know

Hi,

I understand your confusion as this is not clearly explained anywhere in the output. The values within parentheses are those obtained as such from the CPHF/KS calculation. The values outside the parentheses are those obtained with an a posteriori symmetrization.

In your case above, the large components are not affected by the symmetrization, which only helps "cleaning" the vanishing components (i.e. they become cleaner "zeros").

Hope this clarifies things a little,

Hi,

The FINALRUN option works fine in P-CRYSTAL: that sentence in the Manual must be a leftover from the past, sorry for that.

Can you elaborate on how the optimized structures are unreasonable? Is it just for the underestimation of the computed band gap or also for some structural aspects? In this case, it would help if you could show us the "expected" structure, as well as the one you obtain from the optimization.

Hi,

Jefferson Maul has run a few tests and indeed we were able to reproduce the odd behavior you've been experiencing. Those negative values for the electron population of the vacancy are indeed very strange and may be suggestive of a basis set unbalance. It seems that whenever basis functions are left on the vacant site (either with your initial ATOMSUBS approach or with the more "canonical" GHOSTS approach) those functions are much "needed" by surrounding atoms.

We have found a possible way forward for your system through the use of the ATOMREMO option. With this option, you create the vacancy by removing the selected atom alongside its basis functions. For instance, we have used this approach to study Oxygen vacancies in CaSnO3 here.

We have used this approach on your system (on a smaller supercell to be able to run the tests more efficiently) and the geometry optimization went well. See the attached output file as a guide.

Hope this helps,

Hi,

I understand your confusion, which originates from the fact that PRESSURE is an old option of FREQCALC that indeed does not apply any pressure. With this option, the structure is unchanged, the forces are unchanged, and ultimately, the harmonic vibration frequencies are unchanged. The only bit of information that is affected by this keyword is the printed value of the PV term in the thermodynamic analysis.

So, this PRESSURE keyword within FREQCALC should be used in just one way: to set the value of pressure corresponding to the structure used to run the harmonic frequency calculation, so as to have the PV term entering the thermodynamic functions (enthalpy and Gibbs free energy) right. Let me make a couple of examples:

1. You perform a full structural relaxation (atomic positions + lattice parameters). Thus, you have a zero pressure (p=0) structure. Then you compute the harmonic frequencies. If you are interested in thermodynamic potentials (i.e. enthalpy, Gibbs, Helmholtz, etc.) you should use the PRESSURE keyword and set the pressure to zero to get the PV term right in the output (otherwise by default it would be computed at p = 1 atm - do not ask me why).

2. You perform a pressure-constrained geometry optimization with the EXTPRESS keyword of OPTGEOM (let's say at 29.457 GPa as in your example). Then you compute harmonic frequencies with FREQCALC on this compressed structure. Now, again, if you are interested in thermodynamic potentials (i.e. enthalpy, Gibbs, Helmholtz, etc.) you should use the PRESSURE keyword and set the pressure to 29.457 GPa to get the PV term right in the output.

So, as I said at the beginning, the PRESSURE keyword here does not change anything on the structure/forces. It only allows you to tell the program what is the pressure of the structure you are working on, so that the printed PV term is right.

A more formally-sound way to combine temperature and pressure is provided by the quasi-harmonic approximation: See Eq. (13) of the tutorial page on the Quasi-Harmonic Approximation (QHA).

Hope this helps clarifying things a little,

Hi,

When computing (and then plotting) 2D maps of the electron density, or spin density (by use of the ECHG keyword), or electrostatic potential (by use of the POTM keyword), a plane must be specified. Now, in CRYSTAL a plane is specified by three points A, B, and C (i.e. if you specify the coordinates of three -- non rectilinear -- points, you identify a unique plane).

In CRYSTAL, you can specify the three points A, B and C in two ways:

1. You can directly input the Cartesian coordinates of the three points by use of the COORDINA keyword;

2. You can select three atoms as your three points, with the ATOMS keyword. In this case, as the system is periodic and you may need to refer to atoms outside of the reference cell, each atom is identified by four integer numbers (the sequential index of the atom within the reference cell and three indices identifying the cell where your selected atom is centered -- i.e. 4 0 0 1 for instance means atom number 4 in the list, found in cell 0 0 1, that is a cell translated by one lattice vector along the c crystallographic axis).

Now, if the three points A, B and C that you have defined are not "orthogonal", i.e. the angle \( \widehat{ABC} \neq 90^\circ\) this would produce an oblique map. To make it rectangular, you can insert the RECTANGU keyword, as shown below:

You can also change the explored range in the selected plane (relative to that identified by the three original points A, B and C) by use of the keyword MARGINS, as shown below:

Input examples can be found at this tutorial page.

Hope this helps,

Hi,

To run a calculation with a not neutral unit cell, you must insert the keyword CHARGED at the end of the basis set input block, as in:

Title
[geometry]
END
[basis set]
CHARGED
END
[SCF parameters]
END

Hope this helps,

Hi,

You can compute the electron density on a 3D grid of points and obtain data directly saved in a .cube format by use of the ECH3 keyword of the PROPERTIES module. If your system is 3D, all you have to do is insert something like:

ECH3
100
END

in your .d3 PROPERTIES input file, where 100 is the number of points in the grid along the first axis. For a cubic system, this would produce a 100x100x100 grid.

Hope this helps,

Hi,

I haven't used it in a while but I remember XCrySDen had a GUI for k path selection.

What I typically do instead is use the Bilbao Crystallographic Server as a guide. If you navigate to "Space-group Symmetry" and then to "KVEC: The k-vector types and Brillouin zones of Space Groups", you will be able to select your space group of interest and then click on "Comparative listing of k-vector types".

For instance, for your P-1 space group you get:

where you can find the coordinates of the high symmetry special points listed in the table (NOTE: if your space group admits a crystallographic cell that differs from the primitive cell, you will find both sets of coordinates in the table. Make sure you use the primitive ones to write the BAND input in the .d3 file) and also inspect the Brillouin zone:

The order at which you explore these points does not really matter (unless you want to reproduce a specific path from some reference plot), as long as you explore most/all of them.

Hope this helps,

Hi,

Would you please share your full input file?

Dear Aparajita,

Gryffindor said in Problem with restarting CPKS calculation:

I’m reaching out again regarding the CPKS step. After several attempts (and a fair bit of persistence!), I was finally able to complete the SCF calculations.

Good!

Gryffindor said in Problem with restarting CPKS calculation:

As per your previous suggestions, I ensured the SCF was converged beforehand, but unfortunately, the CPKS has never been able to finish successfully on my side.

Is there a reason for switching DIIS off in your CPKS calculation?

Also, you can change the convergence criterion for the CPKS with the TOLALPHA keyword. In your case, it may be enough to run it with:

CPKS
TOLALPHA
2
END

Gryffindor said in Problem with restarting CPKS calculation:

Since I need these results quite urgently, I was wondering if you could try running it on your end to see if there's anything I'm missing?

I am afraid I can not. It is a huge calculation that you are running on 800+ atoms/cell and I do not have the computing power for this.

Dear Aparajita,

It is strange that the fort.9 unit is empty. May it be a problem of the script in copying back that unit from the scratch folder to the working directory?

Anyway, you can try to trick CRYSTAL in this way:

• By inspection of your output file:

CYC  13 ETOT(AU) -1.952535840054E+05 DETOT -1.58E-05

the energy difference is down to 10^-5 at cycle 13. So, you can try to run the calculation with a low energy convergence threshold of just 10^-4 with something like this:

TOLDEE
4

and the SCF would be considered as converged for the program. In this way the calculation is not killed and probably you will get your fort.9 unit,;

• you can then use it to do a restart calculation with the GUESSP keyword (and by tightening the energy convergence threshold to something like TOLDEE 7, or 8):

GUESSP
TOLDEE
7

Hope this helps,

job314 said in No space left on device:

those are huge HPC nodes… they can’t be possibly out of disk…

On our cluster, although the total disk space is huge, there are limited disk quotas for each user. Maybe it is the same there for you?

job314 said in No space left on device:

will it affect my convergence or calculation speed?

Possibly, but not much I think.