Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

Skattejag · 2022-03-19 10:51:56

Dear Prof Tian Lu,

I've been trying to identify and isolate domains from a Gaussian .cube file containing charge density calculation of an inter-fragment/intramolecular interaction, obtained via the molecular tailoring approach. Both positive and negative charge density isosurfaces of the dipolar lobes are visible in Paraview/multiwfn, when appropriate isovalues are set.

When I use subfunction 14 [Domain analysis (Obtaining properties within isosurfaces of a function)] of main function 200, followed by option -1 (Yield domains based on the grid data in memory), after setting the "criterion for defining domain, current: <'to one of the values listed below'>. The calculation seems to run fine for any number of points below 1,000,000 and finishes after a couple of seconds for the Linux code (running on a cluster) and a couple of minutes running on my Windows PC, however when the criteria are set that the number of qualifying grids/points reach above one million, the calculation seems to hang for an undetermined amount of time (see below for several examples).

I do not know if this might be due to the algorithm struggling to separate domains due to many points being to far apart to qualify as being in the same domain and therefore too many domains being found (i.e. my data being the problem)? However, in my opinion, it is suspicious that the problem seems to occur when more than 1,000,000 points are involved. Please let me know if this is a bug or if my data may be to blame.

The following domain analysis searches were run, given in the format: (density) criterion (a.u.) (wall clock time) ... <number of grids/points>:
>0.001 (4 sec) ... 351,037 points or 'grids'
>0.0005 (12 sec) ... 898,559 points or 'grids'

>0.00049 (13 sec) ... 920,337 points or 'grids'
>0.00048 (13 sec) ... 942,620 points or 'grids'
>0.00047 (13 sec) ... 965,990 points or 'grids'
>0.00046 (14 sec) ... 990,287 points or 'grids'
>0.000458 (14 sec) ... 995,368 points or 'grids'
>0.000457 (14 sec) ... 997,877 points or 'grids'

>0.000456 (calculation not finishing after several minutes and in several cases not finishing after several hours) ... 1,000,329 points or 'grids'
>0.000455 (calculation not finishing....) ... 1,002,796 points or 'grids'
>0.00045 (calculation not finishing....) ... 1,015,634 points or 'grids'
>0.0004 (calculation not finishing....) ... 1,160,421 points or 'grids'

Here are the Multifwn (Linux) version details (the Windows version was downloaded at the same time as the Linux version), as well as some details about my input .cube file :

Multiwfn -- A Multifunctional Wavefunction Analyzer
Version 3.8(dev), release date: 2022-Mar-16
Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
Below paper ***MUST BE CITED*** if Multiwfn is utilized in your work:
Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580-592 (2012)
See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information
Multiwfn official website: http://sobereva.com/multiwfn
Multiwfn English forum: http://sobereva.com/wfnbbs
Multiwfn Chinese forum: http://bbs.keinsci.com/wfn

( Number of parallel threads: 24 Current date: 2022-03-19 Time: 07:24:00 )

Input file path, for example E:\Love_Live!\Nico_Yazawa.wfn
(Supported: .mwfn/wfn/wfx/fch/molden/31/chg/pdb/xyz/mol/mol2/cif/cub, etc.)
Hint: Press ENTER button directly can select file in a GUI window. To reload the file last time used, simply input the letter "o". Input such as ?miku.fch can open the miku.fch in the same folder as the file last time used.

Selected file: /mnt/lustre/users/<user>/Multiwfn/Multiwfn_3.8_dev_bin_Linux/2_DPA_41ac_DLPNO_MP2_gas_opt_AVDZ_mdcip_orbopt_A.1_interaction_900.eldens.cube
Please wait...

Title line of this file:
Generated by Multiwfn
Totally 729000000 grid points

Total number of atoms: 28

Grid information:
Translation vector: X Y Z (Bohr)
Vector 1: 0.026126 0.000000 0.000000
Vector 2: 0.000000 0.033761 0.000000
Vector 3: 0.000000 0.000000 0.029878
Volume of each grid: 0.000026 Bohr^3
The range of x is from -9.122529 to 14.364745 Bohr, 900 points
The range of y is from -16.997071 to 13.354068 Bohr, 900 points
The range of z is from -17.764033 to 9.096289 Bohr, 900 points
Total number of grid points: 729000000
This grid data will take up at least 1465 MB memory

Loading grid data, please wait...

Statistical information:
Global minimum value: -0.56700700E+00 at 1.14499 3.32705 -4.67747 Bohr
Global maximum value: 0.64799200E+00 at 1.14499 3.29329 -4.67747 Bohr
Differential element: 0.0000263536 Bohr^3
Summing up positive value in grid file: 3656.4414884829
After multiplied by differential element: 0.0963603515
Summing up negative value in grid file: -3633.3236097234
After multiplied by differential element: -0.0957511124
Summing up all value in grid file: 23.1178787595
After multiplied by differential element: 0.0006092390

Loaded /mnt/lustre/users/<user>/Multiwfn/Multiwfn_3.8_dev_bin_Linux/2_DPA_41ac_DLPNO_MP2_gas_opt_AVDZ_mdcip_orbopt_A.1_interaction_900.eldens.cube successfully!
....

PS: Thank you for developing Multiwfn and for writing such an informative manual. It is such a pleasure to be able to use such a powerful and versatile software.

Thank you for your assistance with this matter.
Kind regards
Cameron

sobereva · 2022-03-19 23:15:10

Dear Cameron,

Usually this problem is caused by inappropriate definition of determining condition of domains. It is best to send me your compressed .cub file via E-mail so that I can check. If the file is too large after compression, please upload to net disk.

Best regards,

Tian

Skattejag · 2022-03-20 08:39:21

Dear Prof Tian Lu,

Thank you for your prompt response.

The Gaussian .cube file is quite large (it's about 10 GB, for 900^3 points). I've made the uncompressed file available via a shared Google drive link (as I don't think compression will be very effective here):
https://drive.google.com/file/d/1ty0eDf … sp=sharing
(Please let me know if you can access this link)

On a separate note:
1) is it possible to export the data for each individual domain as a .pdb or a type of coordinate file (x,y,z,value), instead of only exporting the domain boundary points, in multiwfn. As I would like to 1) integrate the domains, 2) determine the minima/maxima of the charge density and 3) calculate the (atomic?) dipole moments for all the dipole pairs associated with each atom in the system from the Gaussian .cube file.

2) I'm not quite certain how I would go about calculating the aforementioned (atomic?) dipole moments, for the dipoles generated by the molecular tailoring approach.

From section '3.18.3 Calculate atomic and molecular multipole moments and <r2> (2)', in the manual, and from the requirements of an electric dipole moment: 'When two electrical charges, of opposite sign and equal magnitude, are separated by a distance, an electric dipole is established' I have interpreted the information as follows:

The atomic dipole moment is calculated based on the assumption that the 'number of electrons associated with an atom' = 'nuclear charge of an atom' and are defined according to the atomic dipole vector components (involving the relative spatial positions only) that are centred on the nuclear position and terminate at all regions in space, where the atomic dipole vector components may be separated into x,y,z components and mulitplied by the number of electrons - associated with a specific atom- at all points in space, thereafter the sum of all dipole components produce in the x-, y- and z-dimensions may then be converted into the total atomic dipole moment vector, which has a magnitude and direction, by the formula supplied. The molecular dipole moment would then be the sum of all atomic dipole vectors.

For the case of interaction charge given by the molecular tailoring approach, the positive charge would no longer be situated on the nuclear position, but would be delocalised into the + component of the dipole (i.e. the region of electron depletion or negative charge density). How would one compute the dipole moment for this case, perhaps by approximating the dipole moment by determining the vector from the centroid/barycenter of each lobe (from - to + charge density) multiplied by the total number of electrons (minimum?/average?/maximum?) in the lobe(s)?

3) (Please read below for context) Is it possible to generate an 'interaction only' wave function file (eg. molden) (and ultimately obtaining the 'interaction only' electron density data) by applying the below operations to the individual fragment .molden files, i.e. by modifying the MO coefficients, or is it only possible by applying the below operations to the individual fragment Gaussian .cube files? Since the cube files contain Cartesian/spatial data, whereas the .molden files contain atom-centred data.

Here is an example of the molecular tailoring approach for determination of the intramolecular interaction energy/density:
For example, to calculate the intramolecular interaction between two methyl groups in dimethylamine (DMA) in a specific conformation, one would determine the electron density data of
1) dimethylamine (DMA) - of that specific conformation given in Cartessian coordinates- ,
2) the left-hand side(LHS), neutral methylamine (generated by substituting the RHS methyl group with H and performing a cartesian-based constrained geometry optimisation where only the non-physical H is optimised) -> fragment MA1,
3) the right-hand side (RHS) methylamine (generated by substituting the LHS methyl group with H and performing a cartesian-based constrained geometry optimisation where only the non-physical H is optimised)-> fragment MA2,
4) the central NH3 (generated by substituting the both methyl groups with the coordinates of the H atoms as calculated in (2) and (3)) -> fragment NH3

The approximate 'intramolecular charge density' = ρ(DMA) - ρ(MA1) - ρ(MA2) + ρ(NH3), where the scalar field or Gaussian .cube grids - for all components - are determined according to the grid assigned to ρ(DMA).

Kind regards
Cameron

sobereva · 2022-03-20 23:33:15

Dear Cameron,

Please compress your cube file before uploading, the size can be reduced by more than 70%. Also please compress the file as different parts, because in my current network environment I cannot successfully download a very large file from Google drive.

Best regards,

Tian

Skattejag · 2022-03-21 06:21:30

Dear Prof Tian Lu,

Please accept my sincere apologies, here are the files for the compressed Gaussian .cube file split into six 307MB .7z (7-zip) files, giving a total size of 1.74GB [the initial estimation of the compression ratio was 17% (about what I thought), however the final compression ratio was 81.8% (remarkable!)]:
https://drive.google.com/drive/folders/ … sp=sharing

Kind regards
Cameron

sobereva · 2022-03-22 09:58:21

Dear Cameron,

I have updated Multiwfn on official website, now the domain analysis works normally even if number of grids in a domain exceeded 1000000.

In addition, a new option 12 has been added to post-processing menu of domain analysis module, it is able to export X,Y,Z,value of all grids of a selected domain.

I am not sure why you hope to calculate atomic dipole moments based on cube file. They can be directly and analytically calculated based on wavefunction file via option 2 of fuzzy analysis module (main function 15) of Multiwfn.

To calculate 'intramolecular charge density' = ρ(DMA) - ρ(MA1) - ρ(MA2) + ρ(NH3), there is no need to generate an 'interaction only' wave function file. If you have wavefunction files of DMA, MA1, MA2 and NH3, you can directly use the very flexible "custom operation" feature of main function 5 of Multiwfn to calculate grid data of 'intramolecular charge density', which can then also be exported as cube file in post-process menu of main function 5. See Section 4.5.5 of Multiwfn manual for example of using this feature.

Best regards,

Tian

Skattejag · 2022-03-22 19:20:25

Dear Prof Tian Lu,

Thank you so much for updating Multiwfn and for adding option 12, it is much appreciated. It will be much easier for me to extract the required data now.

Sobereva wrote:

I am not sure why you hope to calculate atomic dipole moments based on cube file. They can be directly and analytically calculated based on wavefunction file via option 2 of fuzzy analysis module (main function 15) of Multiwfn.

Please note that I am still familiarising myself with the theory behind determining (atomic) dipole moments. I may not be explaining myself properly.

As I understand it, an atomic dipole moment is determined from the (presence of) electron density data in all regions that are associated with a specific atom's, via discrete or fuzzy partitioning methods, as well as the the nuclear coordinate and nuclear charge. Whereas, the 'intramolecular charge density' consists of regions of space associated with the presence of electron density (+ values/lobes) and the absence of electron density (- values/lobes), relative to the electron density distribution of the conformer/molecule of interest (eg. DMA) and are obtained for a specific set of fragments, i.e. the CH3...CH3 interaction in DMA. For larger molecules (see the .cube file, where >6 sets can be defined), there may be several interfragment sets/pairs that may be defined, where the interaction energy and interaction charge density may be determined for each set. As you can see from the .cube data the dipoles are centred on each atom (astride each atom), but are not spherically distributed around the atoms, so the scenario is different from standard atom-centred GTO data. I'm not certain if these type of interaction-specific atom-centred dipoles are still considered to be a type of atomic dipole.

From the Gaussian .cube file, I suppose I could determine the following contributing components to the (atomic) dipole moment for each interaction set:[+(charge density),+(nuclear)] and [-(charge density),+(nuclear)], where +(nuclear) would be the nuclear charge and where unmerged charge density domains could be analysed for each atom (this would define the atomic space), and then all component dipoles could be combined to generate the total dipole associated with each atom. I don't know if the method/code in the 'Fuzzy atomic space'>'Calculate atomic and molecular multipole moments and <r^2>' module in Multiwfn will give the desired output, these methods utilise generated atomic densities and I do not believe this will model the interaction-specific atom-centred polarised density correctly.

I would like to visualise and evaluate the atom-centred dipole vectors for a chosen interaction for several conformers of a molecule, and I would like to establish what type of correlation exists between the calculated interfragment interaction energies and the 'charge density polarisation'/dipole data.

Sobereva wrote:

To calculate 'intramolecular charge density' = ρ(DMA) - ρ(MA1) - ρ(MA2) + ρ(NH3), there is no need to generate an 'interaction only' wave function file. If you have wavefunction files of DMA, MA1, MA2 and NH3, you can directly use the very flexible "custom operation" feature of main function 5 of Multiwfn to calculate grid data of 'intramolecular charge density', which can then also be exported as cube file in post-process menu of main function 5. See Section 4.5.5 of Multiwfn manual for example of using this feature.

Thank you, I was unaware of the the 'main function 5 > sub-function 0' method of calculating the 'interaction charge density' .cube files. I could generate .molden files, via ORCA, for the molecular fragments and load them using this method. This should speed up the generation of these Guassian .cube files considerably, assuming that 'interaction only' wave function files cannot be generated for analytical evaluation instead.

Please let me know if you have any other questions or suggestions.

Thank you for your assistance with this matter.
Kind regards
Cameron

sobereva · 2022-03-23 02:06:25

Dear Cameron,

Multiwfn calculates atomic dipole moments by integrating product of coordinate vector (w.r.t. atomic nucleus) and electron density within atomic spaces. What you want to obtain, seems to be the variation of atomic dipole moments by integrating product of coordinate vector (w.r.t. atomic nucleus) and 'intramolecular charge density' within atomic spaces, is it correct?

Best regards,

Tian

Skattejag · 2022-03-23 06:31:01

Dear Prof Tian Lu,

sobereva wrote:

Multiwfn calculates atomic dipole moments by integrating product of coordinate vector (w.r.t. atomic nucleus) and electron density within atomic spaces. What you want to obtain, seems to be the variation of atomic dipole moments by integrating product of coordinate vector (w.r.t. atomic nucleus) and 'intramolecular charge density' within atomic spaces, is it correct?

Correct, however I am unsure of how the atomic spaces should be defined without introducing contamination of charge density from other adjacent atoms and ensuring the total associated atomic charge density is fully accounted for. These errors may not be avoidable if a conventional discrete atom partitioning method (eg. AIM basins) or fuzzy atom partitioning methods are used.

Here is an image of six isosurfaces of the interpyridyl-interaction charge density (±0.004, ±0.0004, ±0.00004 a.u.) shown with varying levels of opacity (1.0, 0.3 and 0.1, respectively), generated from the supplied .cube file:
charge bbcode test

From the above, one can see that at high |charge densities| (eg. ±0.004 a.u., the high opacity isosurfaces) the positive- and -negative lobes are clearly delineated and create discrete atom-centred spaces/partitions/domains, however at lower |charge densities| (eg. ±0.00004 a.u., the low opacity isosurfaces) the domains belonging to individual atoms merge. The atomic dipole moments calculated from the latter data would contain contaminated charge density values from other atoms. This is why I proposed determining the atomic dipole moment by specifying an isovalue limit for domain selection (eg. ±0.004 a.u, i.e <-0.004 and >0.004 a.u.).

I hope that I have explained myself clearly and I would also like to thank you for your interest in my problem.

Kind regards
Cameron

sobereva · 2022-03-24 02:45:45

Dear Cameron,

Today I updated Multiwfn. If "ispecial" in settings.ini is set to 1, then the electron density involved in option "2 Calculate atomic and molecular multipole moments and <r^2>" in fuzzy analysis module will be replaced with user-defined function. Via this feature, you can set "iuserfunc" in settings.ini to -1 (user-defined function will correspond to intrapolated function based on grid data, as mentioned in Section 2.7 of manual), then boot up Multiwfn, load the cube file corresponding to 'intramolecular charge density', then use the aforementioned option. Then Multiwfn will calculate atomic dipole/multipole moments based on 'intramolecular charge density'.

I am not sure if this is useful for you, but this what I can do my best currently.

Best regards,

Tian

Skattejag · 2022-03-28 10:42:53

Dear Prof Tian Lu,

My apologies for the delay. Thank you for introducing this new feature in Multiwfn (24/03/2022). This is long post, so feel free to focus your replies to question 1 and 3. The errors shown in point 3 are the most concerning to me at the moment, and I need to decide on which protocol I should use for interaction density .cube file generation :1) Multwfn's 'calc. grid data>set custom operation' route where is specify the .molden files or 2) ORCA's plot generation utility to generate each component .cube file (for each .mdcip file) followed by Multiwfn's 'Grid data. calculation').

1) I have tested the new special feature and it produces output. I would like to know how the atomic space - used in the option "2 Calculate atomic and molecular multipole moments and <r^2>" - and corresponding atomic dipole moments are determined.

Does it use the fuzzy atomic space (i.e. the atomic weighing function) generated via the selected fuzzy partitioning method (for example, the default Becke's method), where it evaluates the trilinearly interpolated .cube data (charge density data) instead of the electron density derived from GTFs - using the aforementioned weighing function- to yield the atomic dipole moments? Does the code take into consideration for charge density contamination originating at atoms other than the selected atom for which the atomic dipole moment is to be determined?

I thank you for your efforts. If I'm not able to use the aforementioned method in my work, I should be able to write a Python script to generate approximate atomic dipole moments via the domain-based method that I proposed, and with the help of the domain export utility that you provided in the domain analysis module. An initial cursory investigation has revealed some potential concerns though, see below.

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

2) For the interpyridyl interaction charge density data: I have evaluated several thresholds of positive and negative domains and I have determined the net charge for each set of domains, via integration of all positive domains and all negative domains, as defined by the domain boundary conditions/thresholds shown in the figure below. It can be seen from the figure that the relative polarisation error [(integral of + density) + (integral of - density )/(integral of + density)*100] contribute significantly for certain domain thresholds. The assumption being that the net charge should be zero for each set of domains, i.e. polarisation of x electrons, should lead to a electron peak of x electrons (+ domain) and a electron hole of x electrons (- domain).

There are two domain thresholds where the relative polarisation errors are ideal or close to 0: -0.002 > ρ > 0.002 and -0.00005 > ρ > 0.00005 a.u. As expected the polarisation error for the full .cube (all data) is close to zero, however significant errors approaching -12.44% and >14.82% can be observed for other domain thresholds. This leads me to conclude that the charge density modelled for the non-ideal domain threshold cases do not show the correct polarisation behaviour. Determination of approximate atomic dipole moments from the | 0.00005 a.u.| data will clearly be problematic, since the atomic charge density is significantly delocalised and merges/mixes significantly with nearby atomic charge density, whereas the atomic dipole moments determined from the |0.002 a.u.| show clearly unmerged/defined atom-centred polarision lobes, which would be ideal for atomic dipole moment determination. The drawback of the latter data being that only a fraction of the total charge density data will be considered, i.e. the |0.002 a.u.| data with a total polarisation of 29.68 milielectrons, compared to the full |0.0 a.u.| data with a total polarisation of 96.06 milielectrons.

charge_density

Perhaps by applying a scaling factor to the sum of the individual atomic dipoles (as determined from the below equations) calculated from the |0.002 a.u.| data to the produce the correct total polarisation for the the full .cube data (i.e. |0.0 (a.u.)| data) can be approximated?

I have written the equations and definitions for the interaction-specific atomic dipole moments that I would like to determine (See below). Please note that I do not have a strong mathematical background, so I may have made several mistakes along the way. From the equations it can be seen that symmetrical polarisation of electron density would produce a norm of the total atomic dipole moment of 0. The atom dipole moment as provided here, would give an indication of the deviation from symmetrical distribution of charge density. Do you have any concerns or suggestions here?

Note: The charge density threshold of + or - 0.00004 a.u. given here is only a placeholder, until I decide how I would like to approach this problem then I will update it.

dipole_eq

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

3) I would like to know if the following discrepancy in electron densities between the generated .cube data given by ORCA and Multiwfn, was due to some setting that I did not adjust correctly in MultiWFN (I did update some values in the settings.ini file, see below for more details)? If not, then was this error due to the conversion of the .optorb file to .molden file format, by the orca_2mkl program? I assume that the ORCA %plots utility correctly determines the electron density that it outputs in the .cube file, when it reads the CCSD MDCI density from the .mdcip file.

I have checked whether the ORCA gaussian .cube files generated match the corresponding Multiwfn-generated .cube files. These data are for the molecular electron densities and not the interaction densities, as I wanted to test/evaluate potential sources of error. I have found the following deviation/error in the number of electrons across all grid points, between the two files (Multwfn .cube relative to ORCA .cube):

(Positive value) After multiplied by differential element: 0.2599997904
(Negative value) After multiplied by differential element: -0.2594284678
After multiplied by differential element: 0.0005713226

The procedure used to generate the two Gaussian .cube files (and the difference .cube) was as follows:

3.a) a .cube was generated, via the %plots module within ORCA 4.2.0, using the CCSD orbitally-optimised .mdcip density.
3.b) a .molden file generated from an orbitally-optimised coupled-cluster wave function file (.mdci.optorb ) was produced after the CCSD calculation finished, via the orca_2mkl program of ORCA 4.2.0., where the .molden file was loaded in Multiwfn (with the following settings.ini options set: 'expcutoff= 1' and 'idelvirorb= 0', in order to generate high accuracy electron density data from the wave-function file). The following sequence of commands were entered after loading the .molden file: '5 Output and plot specific property within a spatial region (calc. grid data)' > '1 Electron density' > '8 Use grid setting of another cube file' > 'here the .cube file generated in (3.a) was supplied' > '2 Export data to a Gaussian-type cube file in current folder'.

A new session of Multiwfn was started, where cube 3.a was loaded and cube 3.b was subtracted from cube 3.a, which was then saved as a new difference .cube file (cube 3.c):
'13 Process grid data' > '11 Grid data calculation' > '4 Subtract a grid file' > 'cube 2.b file path' > '0 Export present grid data to Gaussian-type cube file (.cub)' > 'cube 3.c'.

Next, a new session of Multiwfn was started, where the difference cube 'cube 3.c' was loaded ('2_DPA_41ac_mdcip_orca_min_molden_diff.cube'), which produced the following summary:

Title line of this file:
Generated by Multiwfn
Totally 729000000 grid points
Total number of atoms: 28
Grid information:
Translation vector: X Y Z (Bohr)
Vector 1: 0.026126 0.000000 0.000000
Vector 2: 0.000000 0.033761 0.000000
Vector 3: 0.000000 0.000000 0.029878
Volume of each grid: 0.000026 Bohr^3
The range of x is from -9.122529 to 14.364745 Bohr, 900 points
The range of y is from -16.997071 to 13.354068 Bohr, 900 points
The range of z is from -17.764033 to 9.096289 Bohr, 900 points
Total number of grid points: 729000000
This grid data will take up at least 1465 MB memory
Loading grid data, please wait...
Statistical information:
Global minimum value: -0.59200000E+00 at 1.14499 3.32705 -4.67747 Bohr
Global maximum value: 0.60800000E+00 at 1.14499 3.29329 -4.67747 Bohr
Differential element: 0.0000263536 Bohr^3
Summing up positive value in grid file: 9865.8214309365
After multiplied by differential element: 0.2599997904
Summing up negative value in grid file: -9844.1423103071
After multiplied by differential element: -0.2594284678
Summing up all value in grid file: 21.6791206294
After multiplied by differential element: 0.0005713226
Loaded /mnt/lustre/users/<user>/Multiwfn/Multiwfn_3.8_dev_bin_Linux/2_DPA_41ac_mdcip_orca_min_molden_diff.cube successfully!

I have uploaded the files to a shared Google drive directory (I have compressed all files larger than 300MB, via 7zip):
1) the '2_DPA_41ac_DLPNO_MP2_gas_opt_AVDZ_mdcip_900' file is the .cube generated by ORCA 4.2.0 (using the ELDENSMDCI keyword) (3.a)
2) the '2_DPA_41ac_DLPNO_MP2_gas_opt_CCSD_AVDZ.mdci.optorb_molden' file is the .cube generated by Multiwfn (3.b)
3) the '2_DPA_41ac_mdcip_orca_min_molden_diff' file is the difference cube (3.c)
4) I've also added the .optorb file, the .molden file, .mdcip file and the orca_2mkl binary from ORCA 4.2.0

Here is the link: https://drive.google.com/drive/folders/ … sp=sharing

Thank you for your patience and assistance with this matter.
Kind regards
Cameron

Last edited by Skattejag (2022-04-01 14:35:16)

sobereva · 2022-03-28 11:18:26

Dear Cameron,

Currently I am really too busy to carefully look at questions 2 and 3, my brief answer to 1 and 3:

1: You statement is correct. The "contamination" must occur, as the density in cube file is contributed by all atoms.

3: Last year I conducted a systematic study of wavefunction quality, and I also tested orbital-optimized CCD wavefunction produced by ORCA and used Multiwfn to conduct some calculations based on it, now I share the way I used.

This is an example input file of ORCA

! CCD def2-QZVPP tightscf miniprint
%maxcore  6000
%pal nprocs   18 end
%mdci
density orbopt
NatOrbs true
end
* xyz   0   1
O     -1.20821942   -1.32398795    0.00000000
C     -0.14731692   -0.74598886    0.00000000
C      0.00000000    0.71738053    0.00000000
H      0.80948737   -1.30966046    0.00000000
H     -0.91659640    1.29634172    0.00000000
C      1.20305931    1.28805781    0.00000000
H      2.10557128    0.68624144    0.00000000
H      1.33283869    2.36228399    0.00000000
*

After running, the following script can be used to convert the natural orbitals of corresponding level to .mwfn format, which can be used as input file of Multiwfn for subsequent analysis.

#Convert all .mdci.nat files generated by ORCA CCSD calculation to .mwfn file
#!/bin/bash
icc=0
nfile=`ls *.mdci.nat|wc -l`
for inf in *.mdci.nat
do
((icc++))
echo Converting ${inf} to ${inf//mdci.nat/mwfn} ... \($icc of $nfile\)
mv ${inf} ${inf//mdci.nat/gbw}
orca_2mkl ${inf//.mdci.nat} -molden
Multiwfn ${inf//mdci.nat/molden.input} << EOF > /dev/null
100
2
32
${inf//mdci.nat/mwfn}
0
0
q
EOF
done

You do not need to alter "idelvirorb", as it never affects result of electron density evaluation.

Best,

Tian

Skattejag · 2022-04-04 17:13:03

Dear Prof Tian Lu,

I have been trying to investigate the natural orbital approach that you suggested for determining the orbitally optimised CCD density (No singles excitation contributions are available for the density calculation, so I did not determine OO-CCSD instead I determined OO-CCD), as well as trying to understand the composition of the various ORCA output files (.gbw, .optorb, .nat, .mdcip) in terms of their generated electron density .cube files.

I have received a response from Stanpa (Stan, a PhD Student at Max-Planck-Institut für Kohlenforschung) on the ORCA forums that is a bit concerning (see https://orcaforum.kofo.mpg.de/viewtopic … 16#p38419). Stan says that the only way to obtain the OO-CCD density is via generating the .cube file by using the ELDENSMDCI keyword in the %plots block in ORCA (point 1 below). From your approach it seems as if you have employed point 4 below, which he says will produce the "HF electron density in the optimised OO-CCD basis (orbitals are rotated to make the OO-CCD Lagrangian stationary)." I am quite confused, as there appear to be multiple types of 'HF electron densities'.

So it seems I have to use approach (1) below to generate the OO-CCD electron density .cube files, since I cannot use the generated .molden files (via orca_2mkl).

Here is an excerpt from the Stanpa's response (his responses are in bold, while my questions are in plain text):

stanpa wrote:

skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
1) You are saying that the only way to obtain the OO-CCD electron density information in the format of a Gaussian .cube file is via the .mdcip file (and the ELDENSMDCI keyword, via ORCA)?
Yes
skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
2) Converting the .gbw file to a .molden file (via orca_2mkl), loading it in MultiWFN, and generating a electron density-based Gaussian .cube file will contain the corresponding Hartree-Fock electron density?
Yes
skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
3) Converting the .mdci.optorb file to a .molden file (via orca_2mkl), loading it in MultiWFN, and generating a electron density-based Gaussian .cube file will NOT contain the corresponding OO-CCSD electron density? But will instead contain...what?
It will contain the HF electron density in the optimised OO-CCD basis (orbitals are rotated to make the OO-CCD Lagrangian stationary).
skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
4) Converting the .mdci.nat file to a .molden file (via orca_2mkl), loading it in MultiWFN, and generating a electron density-based Gaussian .cube file will NOT contain the corresponding OO-CCSD electron density? But will instead contain...what?
It will contain the HF electron density in the optimised OO-CCD natural orbital basis (similar to 4, but now the orbitals have been rotated additionally, so that the OO-CCD density matrix is diagonal)
skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
5) Since the CC amplitude information is not available in the .mdci.optorb or .mdci.nat files the electron density cannot be correctly constructed from the orbital information (MO coefficients, GTO information) alone, contained within the aforementioned files?
Correct
stanpa wrote:
The discrepancy comes from the fact that you are comparing electron densities from 2 different levels of theory. The ELDENSMDCI will give you the MDCI density, whereas the .optorb file will give you an SCF level density.
skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
6) So from the above, the Gaussian .cube file generated in scenario (3) would be equal to the Gaussian .cube generated in scenario (2) , even though the .mdci.optorb file was only generated after the OO-CCD calculation reached completion and the .gbw file was generated after the SCF/HF-part completed?
AFAIK, no. The GBW file should contain the orbitals in the canonical basis (Fock matrix is diagonal), which is not the same as the one from (3).
skattejag wrote: ↑Mon Apr 04, 2022 1:24 pm
7) How is the electron density information stored in the binary .mdcip file and how is the electron density value determined for a specific unit volume or voxel at a specific coordinate in 3D space (in the .cube)?
The .mdcip file contains the 1-body density matrix in AO-basis. To understand how cube files are constructed, I would suggest to read some literature.

Thank you for your patience and assistance with this matter.
Kind regards
Cameron

Skattejag · 2022-04-09 21:14:23

Dear Prof Tian Lu,

Continuing from my last post...

I have obtained some Gaussian .cube data from some test calculations with ORCA in order to gain more information about the composition of the .gbw, .mdci.optorb, mdci.nat, .mdcip and .scfp files (please see my new post here : https://orcaforum.kofo.mpg.de/viewtopic.php?f=8&t=8954). From my tests I have found that the .mdci.nat file does indeed produce the same electron density as the .mdcip file (i.e. the OO-CCD electron density), and that the .mdci.optorb file does not produce the HF electron density nor the OO-CCD electron density. I don't know what type of electron density information can be obtained from the .mdci.optorb file.

So your suggestion of using the mdci.nat file for OO-CCD should produce the correct electron density. Unfortunately, I don't have the .mdci.nat files for my big conformers as I didn't specify 'NatOrbs True' in the MDCI block, so I'll be using the Gaussian .cube files generated from the .mdcip files for my work.

Thank you again for your assistance with this matter.
Kind regards
Cameron

sobereva · 2022-04-11 19:19:35

Dear Cameron,

Thank you for your information, I am glad to know that my suggestion is correct and didn't mislead you.

Best,

Tian

Multiwfn forum

#1 2022-03-19 10:51:56

Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#2 2022-03-19 23:15:10

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#3 2022-03-20 08:39:21

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#4 2022-03-20 23:33:15

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#5 2022-03-21 06:21:30

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#6 2022-03-22 09:58:21

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#7 2022-03-22 19:20:25

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#8 2022-03-23 02:06:25

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#9 2022-03-23 06:31:01

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#10 2022-03-24 02:45:45

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#11 2022-03-28 10:42:53

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#12 2022-03-28 11:18:26

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#13 2022-04-04 17:13:03

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#14 2022-04-09 21:14:23

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

#15 2022-04-11 19:19:35

Re: Possible bug in Multiwfn(Win/Linux) 'domain analysis' causing hanging

Board footer