Dear Sayan,

If the system is not flexible, one-stage fitting is adequate.
Even for flexible molecules, if the charges are not used to produce MD trajectory but only used for evaluating couplings based on existing trajectory, two-stage fitting is also not needed.

Yes, chgcons.txt and eqvcons.txt only affect one-stage fitting.

I think Gaussian itself cannot do this, because the code was developed only for single state. To calculate TrEsp, special considerations (e.g. treatment of nuclear charges) should be taken into account.

Best regards,

Tian

Hi Tian,

Thanks. Now it works...!! I was giving the wrong input during the calculation.

Regarding, this two-stage fitting you said for MD simulation purpose it is better. For my work, I want to impose my charges on each atom of a molecule and then calculate the couplings along the trajectory, so should I use two-stage fitting or one stage fitting....?

Also, in the two-stage fitting, I think these chgcons.txt/eqvcons.txt inputs are not taking into consideration.

Another doubt, I have realized that Gaussian itself can calculate the  ESP fitted charges using CHELPG method. So, can Gaussian directly calculate the TrESP or transition Mulliken charges if I requested density=transition=1  pop=(chelpg, dipole)...?  The printed charges at the end are the transition charges or ground state charges...?

Best,
Sayan

Dear Sayan

Regarding "...giving some absurd charges", I am not sure about the reason without more detailed information. In principle, the problem should not be relevant to symmetry of a molecule.

Using the new.molden generated via the N-phenylpyrrole_relaxS2.fch (see final part of Section 4.18.13 of the latest version of manual) as input file, I am able to normally calculate RESP charge using one-stage fitting with default setting, the result is shown below:

 Center Charge 1(C ) -0.226936 2(C ) -0.124563 3(C ) -0.124563 4(C ) -0.226936 5(N ) 0.098841 6(H ) 0.169509 7(H ) 0.139598 8(H ) 0.139598 9(H ) 0.169509 10(C ) 0.333217 11(C ) -0.344827 12(C ) -0.344827 13(C ) -0.037293 14(H ) 0.147376 15(C ) -0.037293 16(H ) 0.147376 17(C ) -0.165887 18(H ) 0.083255 19(H ) 0.083255 20(H ) 0.121593 Sum of charges: 0.000003 RMSE: 0.003014 RRMSE: 0.271635

I am not sure why you failed.

If you need original definition of RESP charge (the main aim is for molecular dynamics simulation), the two-stage fitting should be used. While for other purposes, I suggest only using one-stage fitting, it is more flexible and controllable.

Best,

Tian

Hi Tian,

For my symmetric  molecule, when I am trying to get the ESP fitted charges for the ground state from one stage fitting (option 2) using constraint charge over a few atoms (chgcons.txt), it is giving some absurd charges, however, for the TrESP calculation, the one stage fitting is perfectly fine with chgcons.txt input.

I did a TDDFT calculation, then from the  .fchk file, I generated new.molden (natural orbital using 16 of 200 functions). Then I used this new.moden to get the  ESP fitted charges containing over a few atoms, but it didn't work when I use one stage fitting. However, the  Gaussian cubegen utility works perfectly fine.

Also, always I should have to choose one stage fitting or when we can choose the two-stage fitting....??

Thanks,

Best,
Sayan

Dear sayan,

Because the distribution of fitting grid points is not always in line with point group of symmetric molecule, therefore the atomic charges do not always have symmetry as the molecule itself. To relieve this problem, you can try to increase the number of grid points, or use the eqvcons.txt setting.

I don't quite understand the second question. The transition density is a real space function, clearly it is impossible to be constant for any molecule.

Best,

Tian

Hi Tian,

That was my doubt...!!

However, I am facing a problem regarding the symmetry of a molecule.  I have a perfectly symmetric molecule, but when I am calculating the  ESP fitting charges for the ground state, excited state, and TrESP, the charges are not equivalent for the symmetric atoms.  So, should I have to use eqvcons.txt, in the  RESP module in order to get symmetrical charges distribution....?  If this is the case then why it is not automatically taking the identical symmetrical charge distribution for identical molecules as I am using the optimized geometry. So, every time for symmetrical atoms should I have to take this eqvcons.txt into account as in MD simulation this creates problem...?

Also, is it true if the transition density is constant irrespective of charge distribution for a non-symmetrical molecule, then the excitonic coupling should also be constant...?

Best,
Sayan

Dear Sayan,

Sorry I don't understand your meaning.
Assume that there are N basis functions, then the density matrix in the .fch file essentially is a N*N matrix. Since this is a symmetric matrix, the .fch only records its lower triangle part, namely there are N(N+1)/2 values.

Best,

Tian

Hi Tian,

Thank you very much for your explanation. Now, it's clear to me.

Also, another point which I forget to ask you is how the density matrix field i.e. Total SCF Density corresponding to each basis function in the .fch file is calculated from the symmetrized density matrix...?  I mean how those SCF densities values corresponding to each basis function are obtained from the density matrix...? 

Best,
Sayan

Dear Sayan,

Unlike the transition density matrix, the density matrix of a state is essentially a symmetric matrix, therefore you don't need to consider if symmetrizing it.

If you want to use the ESP fitting charges to approximately evaluate molecular total dipole moment, you should use default setting to derive the atomic charge, namely "Nuclear + Electronic".

However, if what you want is the exact molecular dipole moment directly derived via excited state wavefunction, you should go to subfunction 2 of main function 200, the molecular dipole moment is directly outputted, as shown below:

 Molecular nuclear dipole moment (a.u.): X= 0.000000 Y= -0.000000 Z= -0.000027 Norm= 0.000027 Molecular electron dipole moment (a.u.): X= -0.000000 Y= -0.000000 Z= 0.472688 Norm= 0.472688 Molecular dipole moment (a.u.): X= -0.000000 Y= -0.000000 Z= 0.472661 Norm= 0.472661

Best regards,

Tian

Hi Tian,

Thanks for your clarification. Also one thing, when the natural orbitals are generated from the .fchk file, no symmetrization were done, so symmetrization of density matrix is only required for TrESP charges calculation...right ...?   

During the ESP fitting charges which option I should choose if I want to compare the dipole moment printed at the end for ground state calculation...?

1. Nuclear + Electronic
2. Electronic only

Best,
Sayan

I just updated the version 3.6(dev) manual on Multiwfn website, the Section 4.18.13 now explicitly exemplifies how to perform wavefunction analysis for an excited state based on relaxed density, please have a check.

Dear Sayan,

The Mulliken charges derived based on the NO_0001.molden corresponding to unrelaxed excited state density, however the Mulliken charges directly printed by Gaussian via "TD(Nstates=10, Root=1) Density=Current" corresponds to relaxed excited state density. This point has been described at beginning of Section 4.18.13 using small text size. If you use "density=rhoci" to request Gaussian only generate unrelaxed density, the printed Mulliken charge by Gaussian will be identical to those printed by Multiwfn based on the .molden file.

If you intend to perform wavefunction analysis by Multiwfn for relaxed excited state wavefunction, please notice below information at the beginning of Chapter 4:

Namely you should use "TD(Nstates=10, Root=1) Density" to generate and store relaxed excited state density matrix to .chk file, then convert it to .fch and use it as input file, then use subfunction 16 of main function 200 to convert the density matrix in the .fch to natural orbitals. Finally the resulting .molden file can be used as input file for wavefunction analysis of the 1st excited state.

Best regards,

Tian

Hi Tian,

I generate the natural  orbital for the first excited state and calculate the corresponding Mulliken charges by these steps:

Multiwfn NO_0001.molden
7
5
1

However, the Mulliken charges calculated using Multiwfn for the first excited state and the  Mulliken charges printed in the Gaussian output for input like TD(Nstates=10, Root=1) Density=Current out=wfn, are not the same. My understanding is that if you use TD(Nstates=10, Root=1) Density=Current out=wfn, then the printed Mulliken charges in gaussian output are corresponding to the first excited state, and these charges should be same if I calculate manually by the Multiwfn. However, both I am getting are different.

Regards,
Sayan

Dear Sayan,

1 Yes, you can. In the CHELPG module, you can choose "2 Set grid spacing" and set a value smaller than default to increase the number of fitting points.

2 Please read the Section "4.18.13 Study electronic structure of a single excited state and difference between two excited states" of the latest version of Multiwfn 3.6(dev) manual, this section explicitly describes how to perform wavefunction analysis for specific excited state.

Best,

Tian

Hi Tian,

I need to calculate the ESP and Mulliken charges for the first excited state geometry. So I did a TDDFT calculation by TD(Nstates=10, Root=1) Density=Current out=wfn, then from the .wfn file, I can calculate ESP fitted charges for the first excited state. During this  CHELPG fitting, can I increase the number of more fitting points to get better accuracy...?

Also, I need to calculate Mulliken charges for the first excited state. I have tried with .fchk file but I am getting the Mulliken charges for the ground state using Multiwfn. So, could you please help, how can I get Mulliken charges for the first excited state...??

Kind Regards,
Sayan