How to conduct a stress tensor analysis?

aslozada · 2024-02-10 15:43:45

Hi,
I'm wondering how to conduct a stress tensor analysis within the QTAIM framework using Multiwnf. Any tips?
Checking the Multiwnf manual, I couldn't find anything with those references.

Thanks!

sobereva · 2024-02-10 23:31:24

Its is Multiwfn not Multiwnf.

Currently Multiwfn doesn't formally support stress tensor analysis, you may need to extend source code.

aslozada · 2024-02-11 00:13:12

Thank you very much for your reply.
To confirm, should all extensions via "main function 2" be performed using "userfunc"?

Of course, Multiwfn

sobereva · 2024-02-11 00:58:54

Perhaps you need to write a code to implement a new real space function. You can extend user-defined function by modifying "function userfunc" in function.f90. Then in topology analysis module, you can obtain value of the function at critical points via option 7, or plot it along a topology path.

aslozada · 2024-02-28 21:01:33

I am performing some extension tests to obtain the stress-tensor matrix, see the code (Reference: Using the gradient expansion model. The Journal of Chemical Physics 130, 154104 (2009); doi: 10.1063/1.3098140 ):

!%
!@ ---- Obtain of the electronic stress tensor using
!@ ---- the gradient expansion approximation
!===========================================================
i_matrix = 0.0
st_term = 0.0

i_matrix(1,1) = 1.0; i_matrix(2,2) = 1.0; i_matrix(3,3) = 1.0

! Electron density using a gradient expansion
coeffs(1) = (-2.0 * (3.0 * pi**2)**(2/3))/10.0
coeffs(2) = (-1.0/36.0)
coeffs(3) = (-1.0/18.0)
coeffs(4) = (1.0/12.0)

st_term(:,:,1) = coeffs(1) * elerho**(5/3) * i_matrix
st_term(:,:,2) = 0.0
st_term(:,:,3) = coeffs(3) * laplfac*(elehess(1,1)+elehess(2,2)+elehess(3,3)) * i_matrix
st_term(:,:,4) = coeffs(4) * funchess

!tensor_stress = st_term(:,:,1)+st_term(:,:,2)+st_term(:,:,3)+st_term(:,:,4)
tensor_stress = 0.0 * st_term(:,:,1)+ 0.0 * st_term(:,:,2)+ 0.0 * st_term(:,:,3) + st_term(:,:,4)
!===========================================================
call diagsymat(tensor_stress,eigvecmat,eigval,idiagok) !More robust
if (idiagok/=0) write(*,*) "Note: Diagonization of stress tensor matrix failed!"
write(ifileid,"(' Eigenvalues of stress tensor:',3E18.10)") eigval(1:3)
write(ifileid,*) "Eigenvectors (columns) of stress tensor matrix:"
write(ifileid,"(3E18.10)") ((eigvecmat(i,j),j=1,3),i=1,3)

if (ifuncsel==1) then !Output stiffness of stress-tensor
call sort(eigval)
eigmax=eigval(3)
eigmed=eigval(2)
eigmin=eigval(1)
write(ifileid,"(a,f12.6)") " stiffness of stress-tensor matrix:",abs(eigmin)/abs(eigmax)
write(ifileid,*) "==================================================="

end if

!%
Example: Stress-tensor Stiffness values: BCP in biphenyl model

However, I am interested in implementing the stres-tensor model, describe in The Journal of Chemical Physics 134, 234106 (2011) (see equation above). How could I use any of the second derivative matrices already implemented in Multiwfn for this calculation, especially in the case of r´? I appreciate any suggestion on this issue.

sobereva · 2024-02-28 22:33:50

One-electron density matrix in real space:
γ(r,r')=∑{i}∑{j} χ_i(r) * P_i,j * χ_j(r')
Where i and j loop over all primitive GTFs. P stands for density matrix in primitive GTF basis. χ_i stands for primitive GTF i

To evaluate the equation in your screenshot, you need P (in primitive GTF basis), as well as first and second derivatives of primitive GTFs.

The density matrix in primitive GTF basis can be constructed by "call genPprim", then the global array "Ptot_prim" will be available.

First and second derivatives of primitive GTFs can be easily evaluated, please check "subroutine orbderv". You can find this subroutine loops over all primitive GTFs, and in each loop, the GTFdx, GTFdy, GTFdz evaluated by this subroutine are the three components of first derivative of the present primitive GTF, and GTFdxx, GTFdyy, GTFdzz, GTFdxy, GTFdxz, GTFdyz are components of second derivative. You can properly extract useful piece of code from this subroutine.

aslozada · 2024-03-11 08:56:19

I wrote a patch for the stress tensor matrix calculation by using specific portions of the 'orbderv' subroutine."
The changes are made in the 'function.f90' and 'sub.f90' files.

The following results depict the stiffness and ellipticity calculations of the stress tensor in the biphenyl model. (The wave function has been calculate at the B3LYP/6-311(2d,2p) level)

Within the 'function.f90' module, I added three new subroutines: 'orbderv_stress,' 'gencalcstress,' and 'calc_stress.' After attempting to directly implement the changes in 'orbderv' and 'gencalchessmat,' I opted to reuse portions of those codes in the aforementioned subroutines.

You can find the test patch at this link https://github.com/aslozada/Stress_tens … nsor.patch

Use with Multiwfn_3.8_dev_src_Linux

cp Multiwfn-add-stress-tensor.patch Multiwfn_3.8_dev_src_Linux/

patch -s -p0 < Multiwfn-add-stress-tensor.patch

File to patch: function.f90 ... File to patch: sub.f90

Last edited by aslozada (2024-03-11 21:39:46)

sobereva · 2024-03-11 22:56:38

Great. Please let me know which settings and commands for Multiwfn are needed to reproduce your maps. I will combine your patch into Multiwfn and perform some tests. If the code works well, this analysis will be added to formal release of Multiwfn.

aslozada · 2024-03-12 01:47:38

Dear Prof. Tian,

Thank you very much for your valuable suggestions.

In the current test version of the patch, the calculation is perform for individual coordinates of the CPs using option 7 of the topology module. Plots along the coordinates were generated through a loop in a external script.

Initially, I attempted to utilize the orbderv subroutine directly. However, a persistent memory-related error appeared when adding new arrays in this subroutine. Consequently, I chose to reuse some of the existing code and create a separate subroutine. Indeed, this decision introduces some redundancy in the code, and an improvement needs to be made.

In the current implementation, for the calculation of the stress tensor matrix, direct access to the values GFTdii and GFTdij (i=x,y,z; j=x,y,z) is required.

sobereva · 2024-03-12 23:43:32

OK, when I have spare time recently, I will formally add stress tensor analysis into Multiwfn by adapting your code.

sobereva · 2024-10-20 16:20:51

Stress tensor analysis now is available:
http://sobereva.com/wfnbbs/viewtopic.php?id=1540

Multiwfn forum

#1 2024-02-10 15:43:45

How to conduct a stress tensor analysis?

#2 2024-02-10 23:31:24

Re: How to conduct a stress tensor analysis?

#3 2024-02-11 00:13:12

Re: How to conduct a stress tensor analysis?

#4 2024-02-11 00:58:54

Re: How to conduct a stress tensor analysis?

#5 2024-02-28 21:01:33

Re: How to conduct a stress tensor analysis?

#6 2024-02-28 22:33:50

Re: How to conduct a stress tensor analysis?

#7 2024-03-11 08:56:19

Re: How to conduct a stress tensor analysis?

#8 2024-03-11 22:56:38

Re: How to conduct a stress tensor analysis?

#9 2024-03-12 01:47:38

Re: How to conduct a stress tensor analysis?

#10 2024-03-12 23:43:32

Re: How to conduct a stress tensor analysis?

#11 2024-10-20 16:20:51

Re: How to conduct a stress tensor analysis?

Board footer