MOS in Molden Format

563 views
Skip to first unread message

Luca

unread,
Jul 20, 2018, 4:03:04 AM7/20/18
to cp2k

Dear Developers

I would need to print the Molecular orbitals (MOS) in the molden format  (with cartesian GTOs). 

I checked that the molden format for MOS is available just in the sections MINBAS_ANALYSIS and XAS. 

Could you please give me a few hints on how to modify cp2k for printing the MOS in the molden format?


Best Regards

Luca 

Matt W

unread,
Jul 20, 2018, 4:46:21 AM7/20/18
to cp2k
Hi Luca,

the print key is hidden in the SCF section


it would be good to check it is working really correctly as it hasn't been used much compared to other output formats.

The  code is in molden_utils.F.

Matt

Luca

unread,
Jul 20, 2018, 6:55:45 AM7/20/18
to cp2k
Dear Matt 

Thanks a lot for your answer. Maybe there is a problem in the printing. 
I am not able to understand why some of the coefficients of the GTOs in the Molden Format differ from the ones printed  
in the CP2K format. (please see below) Furthermore the two format differer also on the number of coefficients which are 0.    

FIRST H2O orbital  CP2K - FORMAT:
     1     1  O  2s       -0.830225196
     2     1  O  3s       -0.071821376
     3     1  O  3px       0.000000000
     4     1  O  3py       0.000000000
     5     1  O  3pz      -0.184245524
     6     1  O  4px      -0.000000000
     7     1  O  4py       0.000000000
     8     1  O  4pz      -0.064849590
     9     1  O  3dx2     -0.005672488
    10     1  O  3dxy      0.000000000
    11     1  O  3dxz     -0.000000000
    12     1  O  3dy2      0.003094800
    13     1  O  3dyz      0.000000000
    14     1  O  3dz2      0.002577688

    15     2  H  1s       -0.329191443
    16     2  H  2s       -0.160202370
    17     2  H  2px      -0.000000000
    18     2  H  2py       0.024560726
    19     2  H  2pz      -0.011724829

    20     3  H  1s       -0.329191443
    21     3  H  2s       -0.160202370
    22     3  H  2px       0.000000000
    23     3  H  2py      -0.024560726
    24     3  H  2pz      -0.011724829


First H2O Orbital - MOLDEN FORMAT
 Ene= -0.90227836886466661
 Spin= Alpha
 Occup=   2.0000000000000000
     2  8.302251959088E-01
     3  7.182137568042E-02
     6  1.842455240316E-01
     7  2.417110287291E-12
     9  6.484958964826E-02
    10  5.672487722343E-03
    16  3.291914425302E-01
    17  1.602023703270E-01
    20  1.172482857017E-02
    21  3.291914425209E-01
    22  1.602023703239E-01
    24  2.456072619181E-02
    25  1.172482856855E-02


My input is the following:


&FORCE_EVAL
  METHOD QUICKSTEP
  &DFT
    !WFN_RESTART_FILE_NAME .wfn
    &PRINT
     &MO
      &EACH
       QS_SCF 0
      &END
      FILENAME MOS
      CARTESIAN
      OCCUPATION_NUMBERS
      EIGENVALUES
      EIGENVECTORS
      !MO_INDEX_RANGE 0 5
      NDIGITS 12
     &END
    &END
    BASIS_SET_FILE_NAME cp2k/data/BASIS_SET
    POTENTIAL_FILE_NAME cp2k/data/POTENTIAL
    &QS
      EPS_DEFAULT 1.0E-10
      MAP_CONSISTENT
    &END QS
    &SCF
      MAX_SCF 100
      EPS_SCF 1.0E-6
      SCF_GUESS ATOMIC
      &PRINT
       &MO_ORTHONORMALITY
       &END
       &MOS_MOLDEN
        &EACH
        QS_SCF 0
       &END
      NDIGITS 9
      FILENAME MOS-MOLDEN
     &END
    &END
    &END SCF
    &XC
     &XC_FUNCTIONAL PADE
     &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC 6.0 6.0 6.0
    &END CELL
    &COORD
    O   0.000000    0.000000   -0.065587
    H   0.000000   -0.757136    0.520545
    H   0.000000    0.757136    0.520545
    &END COORD
    &KIND H
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q1
      MAO 5
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q6
      MAO 13
    &END KIND
  &END SUBSYS
&END FORCE_EVAL
&GLOBAL
  PROJECT H2O
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL
~                                                                                                                                          
~                                   

Matt W

unread,
Jul 25, 2018, 6:02:02 AM7/25/18
to cp2k
Hi Luca,

the Molden ordering of p, d, f orbitals is quite different from CP2K.  Details are in the source code or on the Molden webpage.

The NDIGITS keywords gives a cutoff for deciding a coefficient is zero. So if you increase it you should see more entries in the molden output. Setting it small you can get quite small output that is still OK for visualisation.

Matt

Luca

unread,
Jul 28, 2018, 12:05:04 PM7/28/18
to cp2k

Dear Matt

Thank you for your answer. However, I still believe that the printed formats do not differ only for the ordering. Furthermore, modifying NDIGITS does not fix the issue. 

This means that the printed MOLDEN file cannot be used for quantitative post-processing.

I performed tests on H atoms and the H2 molecule, and in both the cases, the MOLDEN-file does not contain the right orbitals. 

Below I will report the H case. The orbital printed in the MOLDEN FORMAT is not normalized, as you can easily check by printing the OVERLAP matrix. 

Please find below the different formats and the overlap (NDIGITS = 12) :

 

CP2K-FORMAT : Occupied orbital of H

     1     1  H  1s        0.860523926375    

     2     1  H  2s       -0.147897200557    

     3     1  H  2px      -0.000000000000   

     4     1  H  2py      -0.000000000000   

     5     1  H  2pz      -0.000000000000  

 

MOLDEN-Format :Occupied orbital of H

Ene= -0.25103226963873182

 Spin= Alpha

 Occup=   1.0000000000000000

     2  8.605239264E-01

 

OVERLAP

                                        1s                               2s 

1    1  H     1s        1.000000008914    -0.933553473764

2    1  H     2s       -0.933553473764     1.000000029566

-----------------------------------------------------------------------------------------------------------------------------------------

 

The right  MOLDEN FORMAT should be (I tested it):

     1    0.860523926375    

     2   -0.147897200557   

 

Best Regards.

Luca 

Matt W

unread,
Jul 28, 2018, 12:41:21 PM7/28/18
to cp2k
Hmmm, could you try recompiling after changing the lines around 323 - 330 to something like 

     DO orbital = 1, 15
         IF (orbmap(orbital) .NE. 0) THEN
            WRITE (iw, fmtstr1) irow_in+orbital, mo_coeff(orbmap(orbital))
            IF (mo_coeff(orbmap(orbital)) .GT. 10.0**(-digits)) THEN
                WRITE(6,*) "skipping this orbital"
            END IF
         END IF
      END DO

i.e. take the original write statement out of the if clause. This should print all coefficients in molden format. 

If that works then I can try and see why the if statement doesn't work as expected ...

Matt

On Saturday, July 28, 2018 at 5:05:04 PM UTC+1, Luca wrote:

Dear Matt

Thank you for your answer. However, I still believe that the printed formats do not differ only for the ordering. Furthermore, modifying NDIGITS does not fix the issue. 

This means that the printed MOLDEN file cannot be used for quantitative post-processing.

I performed tests on H atoms and the H2 molecule, and in both the cases, the MOLDEN-file does not contain the right orbitals. 

Below I will report the H case. The orbital printed in the MOLDEN FORMAT is not normalized, as you can easily check by printing the OVERLAP matrix. 

Please find below the different formats and the overlap (NDIGITS = 12) :

 

CP2K-FORMAT : Occupied orbital of H

     1 1 H 1s 0.860523926375    

     2     1  H  2s       -0.147897200557    

     3 1 H 2px -0.000000000000   

     4     1  H  2py      -0.000000000000   

     5     1  H  2pz      -0.000000000000  

 

MOLDEN-Format :Occupied orbital of H

One = -0.25103226963873182

Luca

unread,
Jul 30, 2018, 2:30:38 AM7/30/18
to cp2k

Dear Matt


The following code seems to give the correct MOLDEN-FORMAT.  Please note that: 

 1) I modified  (irow_in+orbital ) in  (irow_in+orbital – 1)

 2) I tested molden formatted files containing no more than s,p,d functions.

 

DO orbital = 1, 15

   IF (orbmap(orbital) .NE. 0) THEN

       WRITE (iw, fmtstr1) irow_in+orbital - 1 , mo_coeff(orbmap(orbital))

   END IF

END DO

 

By adding the following clause:

  IF (mo_coeff(orbmap(orbital)) .GT. 10.0**(-digits))

 

One obtains an incomplete MODLEN formatted file since it misses some coefficients greater than  10.0**(-digits).

 

Best Regards.

Luca 

Luca

unread,
Jul 31, 2018, 2:09:31 AM7/31/18
to cp2k
Dear Matt

the correct clause should be

 IF (ABS(mo_coeff(orbmap(orbital))) .GT. 10.0**(-digits))


Best Regards
Luca

Luca

unread,
Aug 4, 2018, 5:58:17 AM8/4/18
to cp2k

Dear Matt


Could you please let me know if the orbitals projected on the CARTESIAN - GTOs are normalized?


The basis set (Normalized-Cartesian-GTOs) printed in the MOLDEN formatted file is consistent with the basis set printed by the NWCHEM code. 

When I use the MOLDEN formatted file printed by NWCHEM, my code gives the right number of electrons. 

When I use the orbitals printed by CP2K, I am not able to get the right number of electrons in the system. 

Considering that the basis set is the same, the CP2K orbitals projected on the Normalized-Cartesian-GTOs could not be normalized.


Best Regards
Luca 

Matthias Krack

unread,
Apr 22, 2021, 3:33:45 PM4/22/21
to cp2k
Dear Luca

I have just applied a fix proposed by Hossam Elgabarty to the CP2K trunk (master) affecting the transformation of spherical to Cartesian orbitals which includes also the MOLDEN output. Perhaps that fix resolves the problem which you observed.

Best regards

Matthias

Luca

unread,
Apr 26, 2021, 4:18:43 AM4/26/21
to cp2k
Dear Matthias

thank you a lot for the news.

Best Regards 
Luca 

Luca

unread,
Sep 1, 2021, 4:02:39 PM9/1/21
to cp2k
Dear Matthias

Which line do I have to modify for printing the basis set (exponents and contraction coefficients) in the molden file with more than 6 digits?

Best Regards
Luca 

Krack Matthias (PSI)

unread,
Sep 2, 2021, 4:51:11 AM9/2/21
to cp...@googlegroups.com

Dear Luca

 

Change the format in line 157 of cp2k/src/molden_utils.F from FMT="((T51,2F15.6))" to FMT="((T51,2F20.12))" for instance. If that works fine with MOLDEN then we could apply the NDIGITS keyword also to that format.

 

Best regards

 

Matthias

 

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/0c10e020-799a-4c83-95aa-9fe8a7d64c86n%40googlegroups.com.

Luca

unread,
Sep 2, 2021, 7:16:58 AM9/2/21
to cp2k
Dear Matthias
Yes, changing that line changes the number of basis set digits in the molden file in the right way. 

Thank you.
Best Regards.
Luca 

Krack Matthias (PSI)

unread,
Sep 2, 2021, 10:32:12 AM9/2/21
to cp...@googlegroups.com

Dear Luca

 

FYI, in the CP2K trunk version the keyword NDIGITS controls now also the accuracy of the MOLDEN basis set output.

 

Best regards

 

Matthias

 

Reply all
Reply to author
Forward
0 new messages