[QE-users] Format of the interatomic force constant .fc file

Mon May 12 09:11:53 CEST 2025

Hello,

I have written a subroutine that computes the forces from force 
constants as a function of atomic displacements from their equilibrium 
positions as part of a TDEP implementation here: 
https://github.com/anharmonic/d3q/blob/main/thermal2/harmonic.f90

It is a bit more complicated because it computes the force on atoms in a 
supercell n×m×l starting from force constants for the unit cell from a 
q-point grid n×m×l, which requires a bit of mapping to decide which atom 
is which. I.e. the force acting on atom i in cell R along direction x is

f_ixR = - \sum_jyR' F(ix0; jyR'-R) u_ijR'

Where u_jyR' is the displacement in Bohr f the atom y in cell R' from it 
equilibrium position along direction y. One has to take into account 
that the force constants are only stored with R=0, so on has to use 
translational invariance to get the other elements.

You can find more information in Appendix A of this manuscript: 
https://people.impmc.fr/lpaulatto/HDR.pdf

You could also probably find something similar, but in python, as part 
of the SSCHA package https://sscha.eu/ (or it could also be in Fortran, 
I'm not sure what they have re-implemented when moving to python)

hth

On 11/05/2025 19:29, Chirantan Pramanik wrote:
> Dear Dr. Lorenzo,
>
> Thank you for your reply and the code. I am not familiar with 
> Fortran. I need to calculate the total force constant of the atom 
> numbered 205 (Mg). I can understand that each pair of atoms has a 
> force constant matrix with 9 elements, and for my case, for one 
> supercell, there are a total of {(205*2)*9}=3690 values of matrix 
> elements only related to the 205^th   atom. I need the addition of all 
> of these from each pair of root mean square values. Is there any 
> direct way or any calculation in QE itself to find this?
>
> Thanks,
> Chirantan
>
> Chirantan Pramanik
> Postdoctoral Researcher
> Dept. of Earth and Planetary Sciences
> Weizmann Institute of Science
> Rehovot, Israel
> ------------------------------------------------------------------------
> *From:* users <users-bounces at lists.quantum-espresso.org> on behalf of 
> Lorenzo Paulatto <lorenzo.paulatto at cnrs.fr>
> *Sent:* Tuesday, May 6, 2025 11:09 AM
> *To:* users at lists.quantum-espresso.org <users at lists.quantum-espresso.org>
> *Subject:* Re: [QE-users] Format of the interatomic force constant .fc 
> file
> *Caution: External Sender. Do not click on links or open attachments 
> unless you recognize the sender.*
>
>
> On 05/05/2025 12:52, Chirantan Pramanik wrote:
>> Dear All,
>>
>> I am stuck with the format of the interatomic force constant file 
>> obtained by the q2r.x code. Below is a preview of the .fc file. 
>> Please help me to understand the entries. I need the interatomic 
>> force constant, suppose, only for the atom '205'.
>
> Hello,
>
> please note that force constants always act between two atoms, 
> assuming you know what you are doing, here is the subroutine that 
> reads the force constants stripped of all the parallel distribution 
> and ssafety checks, it is only a couple dozen liens of code. The 
> subroutine reads the force constants for a supercell of size 
> nr1×nr2×nr3, because of translational invariance, the first atom 
> (i,na) is always in the first unit cell (i.e. that with coordinates 
> 0,0,0) while the second atom (j,nb) is in cell (nr1-1)*a1 + (nr2-1)*a2 
> + (nr3-1)*a3. (Or maybe it is the second atom which is in cell 0,0,0, 
> I never remember).
>
> Because in your case the cell is 1×1×1, it won't make any difference.
>
>
> hth
>
>
>
>
> !-----------------------------------------------------------------------
> SUBROUTINEreadfc( flfrc, nr1, nr2, nr3, epsil, nat, &
> ibrav, alat, at, ntyp, amass, omega, &
> has_zstar, alph, read_lr)
> !-----------------------------------------------------------------------
> !
> IMPLICITNONE
> ! I/O variable
> CHARACTER(LEN=256) ::flfrc
> CHARACTER(LEN=80) ::line
> INTEGER::ibrav, nr1,nr2,nr3,nat, ntyp
> REAL(DP) ::alat, at(3,3), epsil(3,3), alph
> LOGICAL::has_zstar, read_lr
> ! local variables
> INTEGER::i, j, na, nb, m1,m2,m3
> INTEGER::ios, ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid
> REAL(DP) ::amass(ntyp), amass_from_file, omega
> INTEGER::nt
> !
> OPEN(unit=1,file=flfrc,status='old',form='formatted')
> !
> ! Read cell parameters
> READ(1,*) ntyp,nat,ibrav,(celldm(i),i=1,6)
> if(ibrav==0) then
> read(1,*) ((at(i,j),i=1,3),j=1,3)
> end if
> !
> CALLlatgen(ibrav,celldm,at(1,1),at(1,2),at(1,3),omega)
> alat =celldm(1)
> at =at / alat ! bring at in units of alat
> CALLvolume(alat,at(1,1),at(1,2),at(1,3),omega)
> !
> ! read atomic types, positions and masses
> !
> ALLOCATE(atm(ntyp))
> DOnt =1,ntyp
> READ(1,*) i,atm(nt),amass_from_file ! type and mass
> END DO
> !
> ALLOCATE(tau(3,nat), ityp(nat), zeu(3,3,nat))
> !
> DOna=1,nat
> READ(1,*) i,ityp(na),(tau(j,na),j=1,3) ! positions of atoms in the cell
> END DO
> !
> ! read macroscopic variables
> ! (dielectric constant and effective charges)
> !
> alph =1.0_dp
> READ(1,'(a)') line
> READ(line,*,iostat=ios) has_zstar, alph
> !
> IF(has_zstar) THEN
> READ(1,*) ((epsil(i,j),j=1,3),i=1,3)
> DOna=1,nat
> READ(1,*)
> READ(1,*) ((zeu(i,j,na),j=1,3),i=1,3)
> END DO
> ELSE
> zeu (:,:,:) =0.d0
> epsil(:,:) =0.d0
> END IF
> !
> READ(1,*) nr1,nr2,nr3
> !
> ! read real-space interatomic force constants
> !
> ALLOCATE( frc(nr1,nr2,nr3,3,3,nat,nat) )
> frc(:,:,:,:,:,:,:) =0.d0
> ALLOCATE( frc_lr(nr1,nr2,nr3,3,3,nat,nat) )
> frc_lr(:,:,:,:,:,:,:) =0.d0
> DOi=1,3
> DOj=1,3
> DOna=1,nat
> DOnb=1,nat
> READ(1,*) ibid, jbid, nabid, nbbid
> READ(1,*) (((m1bid, m2bid, m3bid, &
> frc(m1,m2,m3,i,j,na,nb), &
> m1=1,nr1),m2=1,nr2),m3=1,nr3)
> END IF
> END DO
> END DO
> END DO
> END DO
> !
> CLOSE(unit=1)
> !
> RETURN
> END SUBROUTINEreadfc
>
>
>
> -- 
> Dr. Lorenzo Paulatto
> IdR @ IMPMC - CNRS UMR 7590 & Sorbonne Université
> phone: +33 (0)1 442 79822 / skype: paulatz
> http://people.impmc.fr/lpaulatto/ <http://people.impmc.fr/lpaulatto/> 
> - https://anharmonic.github.io/ <https://anharmonic.github.io/>
> 23-24/423 B115, 4 place Jussieu 75252 Paris CX 05
>
> _______________________________________________________________________________
> The Quantum ESPRESSO Foundation stands in solidarity with all civilians worldwide who are victims of terrorism, military aggression, and indiscriminate warfare.
> --------------------------------------------------------------------------------
> Quantum ESPRESSO is supported by MaX (www.max-centre.eu)
> users mailing listusers at lists.quantum-espresso.org
> https://lists.quantum-espresso.org/mailman/listinfo/users
-- 
Dr. Lorenzo Paulatto
IdR @ IMPMC - CNRS UMR 7590 & Sorbonne Université
phone: +33 (0)1 442 79822 / telegram: lpaulatto
http://www.impmc.upmc.fr/~paulatto/ - https://anharmonic.github.io/
23-24/423 B115, 4 place Jussieu 75252 Paris CX 05
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20250512/0d673a15/attachment.html>