Multiwfn forum

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.

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.