Hi,

I am doing a frequency with intensities calculation on an Iron complex with four Fe atoms in the unit cell. I am trying to do the ferromagnetic state with HSE06-D3 functional and the POB-tzvp-rev2 basis. I have optimised the cell, copied the optimised cell to fort.34 and copied the fort.9 file of the optimised calculation to fort.20, ready for a frequency calculation with intensities.

The calculation is being run on 24 cores on a single node, with up to 4GBs of memory /core. The run uses intel mpi for interprocess communication.

After reporting the hyperpolarisability and electric susceptibility tensors (Beta and Chi), the program crashes with:

COMPUTING IR TENSOR ALONG DIRECTION X forrtl: severe (67): input statement requires too much data, unit 2, file /mnt/scratch/chmjke/SpinCrossover/CRYSTAL/HSE06/AFM2/fort.2.pe15

I attach the input file: frequencies.d12
A tar file of the output and fort.34: info.tgz

Any help would be grateful received.
Thanks John

Let me just add that we do have a development version of the code for HSE+SOC, which we plan to include in the next release.

Hello!
I found the source file. I have also learned that the .f9 was calculated with Crystal17, and that malloc happened in Props17, not in Props23 (I am likewise not sure whether the calculation that failed for my colleague was done using v17 or v23). This might have been a limitation of the machine that I was using, because the same calculation didn't crash with NEWK 36 72.
malloc.d12

I am sorry if that turns out irrelevant, and of no use in the current version of Crystal.

Hi,

Let me just add that of course in P1 there are no special high-symmetry points to guide you in the definition of the path, so you need to be a little creative. For instance, you can start from Gamma (0 0 0) and go to the edge of the FBZ along the b1 reciprocal lattice (1/2 0 0), to then go to (1/2 1/2 0), then to (0 1/2 0) then back to Gamma (0 0 0) and then to the edge along the b3 reciprocal lattice (0 0 1/2). Or something else! 🙂

Hi andrejsc,
There is no "hardcoded" switch in the TOLINTEG keyword, I suspect that the quote from the manual comes from the "old days", when a threshold of \( 10^{-20} \) was considered absurdly small. However, with advances in hardware and the increasingly complex structures we want to compute, such thresholds are sometimes necessary.

A small note: don't warry about going past machine precision with these small numbers. All evaluations of integral thresholds in the code are performed at the logaritmic level (ie only on the exponent), so you should be fine even with a value of 1 million in TOLINTEG (though maybe not fine in terms of computation time)

If you fear that a reduced symmetry might influence the thermodynamic results, you can simply run the calculation without imposing any symmetry at all. For a 3-atom system like CO₂ the computational cost is so minimal that you can run both symmetry-free and \(D_{4h}\) calculations. Comparing the two results should give you an indication of whether the imposed symmetry has any impact on the quantities you are interested in.

Hi Orion,
Looking at your input, I think you want to plot the orbital projected electronic DOS (in CRYSTAL this is done with the DOSS keyword, PDOS in CRYSTAL normally refers to phonon DOS).

To obtain separate px py pz projections, the correct approach is:
set the first parameter of DOSS to 2 (because you want two separate projections).
After the main DOSS line, you must then add one line per projection, specifying which atomic orbitals (AOs) belong to each projection.
Each line has the form:

n AOs index_1 index_2 ... index_n

The AO indices correspond to the ones printed in the CRYSTAL output. To find them, search in the SCF output for "LOCAL ATOMIC FUNCTIONS BASIS SET", you will see a table with all the basis set, that will look like this:

******************************************************************************* ATOM X(AU) Y(AU) Z(AU) N. TYPE EXPONENT S COEF P COEF D/F/G COEF ******************************************************************************* 1 O -2.793 -4.838 -4.110 1 S 8.589E+03 1.895E-03 0.000E+00 0.000E+00 1.297E+03 1.439E-02 0.000E+00 0.000E+00 2.993E+02 7.073E-02 0.000E+00 0.000E+00 8.738E+01 2.400E-01 0.000E+00 0.000E+00 2.568E+01 5.948E-01 0.000E+00 0.000E+00 3.740E+00 2.808E-01 0.000E+00 0.000E+00 2- 5 SP 4.212E+01 1.139E-01 3.651E-02 0.000E+00 9.628E+00 9.208E-01 2.372E-01 0.000E+00 2.853E+00-3.274E-03 8.197E-01 0.000E+00 6- 9 SP 9.057E-01 1.000E+00 1.000E+00 0.000E+00 10- 13 SP 2.556E-01 1.000E+00 1.000E+00 0.000E+00 14- 18 D 1.292E+00 0.000E+00 0.000E+00 1.000E+00

The indeces are the numbers befor each orbital shell (1 2-5 6-9 10-13 14-18).
The notation 2-5 means indices 2 through 5, because that shell contains 4 AOs.
From these blocks you must identify the p-type orbitals for each atom (Zn and Cl). The p shells will appear in groups of three basis functions (px, py, pz).

So, in your system:

Locate the Zn atom in this printed basis list Identify the block corresponding to its p orbitals Extract the AO indices for px, py, pz Do the same for Cl

Then you define one projection per line. For example, structurally something like:

NEWK 12 12 1 0 DOSS 2 100 3 6 1 12 0 3 <list of Zn p AOs> 3 <list of Cl p AOs> END

Keep in mind that, if your basis is double- or triple-Z, each p shell appears as a separate triplet of AOs (one triplet per Z). To correctly project all p you must include the corresponding AOs from every p-shell triplet for that atom.

Let me know if you manage to do this, or if you need further help.

Dear Jonas,

I tried using pymatgen to extract the point-symmetry information from your .xyz file (see the Python script below):

from pymatgen.core import Molecule from pymatgen.symmetry import analyzer bigstructure = Molecule.from_file("yourfile.xyz") PGstructure = analyzer.PointGroupAnalyzer(bigstructure) sym_mol = PGstructure.get_equivalent_atoms() print(sym_mol["eq_sets"])

This returns a Python data structure containing the symmetry-irreducible sets of atoms (only 6 for this system, which is indeed of Ih point group!).
When preparing the CRYSTAL input, be careful with the orientation of your asymmetric unit. In my case, for example, I had to change the sign of the x and y coordinates to make the symmetry consistent with CRYSTAL’s conventions.
Icosahedral point groups are available in CRYSTAL (Ih is point group number 47 in CRYSTAL), so the input fort this molecular cage reads:

Symm. structure MOLECULE 47 6 8 2.605032231 -11.914271806 11.762689798 6 4.344538598 15.664236912 4.366798884 6 3.427479683 14.428996321 8.326818862 6 -8.906580884 1.370529673 14.411150640 1 2.632960331 14.962480562 7.832453804 5 -8.776255231 -2.916256114 14.200279302 COORPRT TESTGEOM END

aerba Christmas is already in the air indeed!

Alright, it sounds like the difficulties originate from the system itself, rather than the particular CRYSTAL input.

Thanks again for all the help!

Best regards,
Wim

aerba, Thank you very much for your reply!

Unfortunately, I don't speak English very well, which is why I may have phrased my question incorrectly.

My question was that the function OPTGEOM does not work correctly with RUNCONFS. I had already tried the steps you suggested.

However, I was able to figure out the problem. For geometric optimization, I just needed to specify the keyword more correctly:

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 OPTGEOM MAXCYCLE 500 ENDOPT END END

I'll leave it here, maybe it will be useful to someone in the future.

Nevertheless, thank you very much for your responsiveness

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,

Good morning, OH that’s cool ☺️, here we go and that is exactly what I want…

Thank you again Dr. Casassa

Thanks. SMEAR helped me to converge scf.

Dear George,
According to your output, you are using CRYSTAL17 to perform the HFsol-3c calculation.
Please note that HFsol-3c and the other "sol-3c" composite methods are only available in CRYSTAL23.

Best regards,
Bartolomeo

Hi Prof. Erba,

Thanks for the reply. From the file, I could see there are value from energy " -7.9324E-01" (Line No: 81826). I could plot the data using Python, Xmgrace. With CRYSPLOT somehow the data loading is slow and my browser hangs so I am not able to check it. But definitely there is data.

I requested for COHP between Ni and three nearby oxygen. In the file there are 5 columns, so if I am not mistaken it corresponds to "Energy - Next three columns in the order I request - Final Column of DOS ?" (Because the last column data when plotted matches the DOS)

Thanks & Regards,
Rams

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,

Thank you, Dr. Alessandro Erba, for the helpful suggestion regarding the use of the COORDINA option. Your example for defining the XY plane was clear and easy to follow — I’ll definitely apply this approach in my setup.

Thanks again for your guidance!
Cheers

Everything works with mpiifx, thank you very much!

Hi Giacomo,

Many thanks for your answer and suggestions.
It works and will do a double check for the geometry.

Have a good day,
Ilias

Hi Chris,

The reason you only see a few k-points printed is indeed that the printing options are not fully supported in parallel execution. In a parallel run, each MPI process handles a subset of the k-points, and only process 0 writes to the output file. As a result, you’ll only see the k-points assigned to process 0, which can be just a few, or even none, depending on how they are distributed.

By default, the code distributes k-points in a round-robin fashion starting from the last process, so the gamma point is typically assigned to the last process and therefore never printed in a parallel execution.

This behavior applies to both options 66 and 67.

Try running the following input:

NEWK 4 4 1 2 66 999 67 999 END

However, you’ll need to run it on a single process for the eigenvectors to be printed correctly.

If the calculation is too expensive to run on a single process, we may need to find an alternative approach to extract those values (though that won’t be straightforward).