Issue with extracting 4th-order force constants

[ZhunyunTang]

Dear Developers,

While using extract_forceconstants, I have encountered a puzzling technical issue, and I kindly ask for your advice and help.

Problem description:
I am trying to extract 4th‑order force constants from 100 MD trajectory snapshots of a supercell containing 375 atoms.

When I set the cutoff radius to 3.5 Å, the program finishes quickly and successfully, producing a force‑constant file of only about 1 MB. This suggests that the actual number of high‑order force constants that satisfy symmetry and distance constraints is not huge.

However, when I increase the cutoff radius only slightly to 4.0 Å, the program crashes ：
(Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x14a1b9e475af in ???
#1 0x5596c427390e in ???
#2 0x5596c428df7f in ???
#3 0x5596c3fbc84a in ???
#4 0x5596c3e6aab9 in ???
#5 0x5596c3e77f72 in ???
#6 0x14a1b9e337e4 in ???
#7 0x5596c3e64b58 in ???
#8 0xffffffffffffffff in ???
./run.sh: line 4: 3872592 Segmentation fault (core dumped) extract_forceconstants -rc4 4.0 -s 1 > fc4.log).

Notably, I ran the test on a large‑memory node with 1 TB of RAM, and I have monitored memory usage in real‑time. No out‑of‑memory condition was observed. Does this mean the problem might stem from internal array allocation strategies rather than the final storage size.
Given that my computational resource is not modest (a 1 TB memory node), the failure caused by a small increase in cutoff radius is very puzzling.
Thank you again for your time.
fc4.log

Best regards,
Zhunyun

[mjv500]

Hello Zhunyun, I suspect your analysis is correct, it's not a memory issue. Could you 1) recompile with debug options, so at least "-g -Wall -fbounds-check -ffpe-trap=invalid,zero,overflow" Then rerun and see if it gives more useful output 2) run under gdb:

gdb `which extract_forceconstants`

gdb> run -rc4 4.0 -s 1 ... and see where it crashes. It should give you even more information. This will pinpoint whether it is doing an allocation, calling an array etc... Best Matthieu Verstraete

…

On Fri, May 29, 2026 at 7:11 AM Zhunyun Tang ***@***.***> wrote: *ZhunyunTang* created an issue (tdep-developers/tdep#178) <#178> Dear Developers, While using extract_forceconstants, I have encountered a puzzling technical issue, and I kindly ask for your advice and help. Problem description: I am trying to extract 4th‑order force constants from 100 MD trajectory snapshots of a supercell containing 375 atoms. When I set the cutoff radius to 3.5 Å, the program finishes quickly and successfully, producing a force‑constant file of only about 1 MB. This suggests that the actual number of high‑order force constants that satisfy symmetry and distance constraints is not huge. However, when I increase the cutoff radius only slightly to 4.0 Å, the program immediately crashes ： (Program received signal SIGSEGV: Segmentation fault - invalid memory reference. Backtrace for this error: #0 0x14a1b9e475af in ??? #1 <#1> 0x5596c427390e in ??? #2 <#2> 0x5596c428df7f in ??? #3 <#3> 0x5596c3fbc84a in ??? #4 <#4> 0x5596c3e6aab9 in ??? #5 <#5> 0x5596c3e77f72 in ??? #6 <#6> 0x14a1b9e337e4 in ??? #7 <#7> 0x5596c3e64b58 in ??? #8 <#8> 0xffffffffffffffff in ??? ./run.sh: line 4: 3872592 Segmentation fault (core dumped) extract_forceconstants -rc4 4.0 -s 1 > fc4.log). Notably, I ran the test on a large‑memory node with 1 TB of RAM, and I have monitored memory usage in real‑time. No out‑of‑memory condition was observed. Does this mean the problem might stem from internal array allocation strategies rather than the final storage size. Given that my computational resource is not modest (a 1 TB memory node), the failure caused by a small increase in cutoff radius is very puzzling. Thank you again for your time. fc4.log <https://github.com/user-attachments/files/28377171/fc4.log> Best regards, Zhunyun — Reply to this email directly, view it on GitHub <#178?email_source=notifications&email_token=ABWOT33WQNHBEDVCLMFKTI345ELWVA5CNFSL4Z3JMQ5C6L3HNF2C22DVMIXUS43TOVSS6NBVGQ2TQNBRG4YDFJTSMVQXG33OVJZXKYTTMNZGSYTFMSSWK5TFNZ2KYZTPN52GK4S7MNWGSY3L>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABWOT3ZXUXEGNNBNQTWND3D45ELWVAVCNFSM6AAAAACZR4TBI6VHI2DSMVQWIX3LMV43ASLTON2WKOZUGU2DKOBUGE3TAMQ> . Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS <https://github.com/notifications/mobile/ios/ABWOT36RFAMRKYK6HMBYGOD45ELWVA5CNFSL4Z3JMQ5C6L3HNF2C22DVMIXUS43TOVSS6NBVGQ2TQNBRG4YDFJTSMVQXG33OVJZXKYTTMNZGSYTFMSSWK5TFNZ2KUZTPN52GK4S7NFXXG> and Android <https://github.com/notifications/mobile/android/ABWOT3Z5AENBVRCBPKW5BI345ELWVA5CNFSL4Z3JMQ5C6L3HNF2C22DVMIXUS43TOVSS6NBVGQ2TQNBRG4YDFJTSMVQXG33OVJZXKYTTMNZGSYTFMSSWK5TFNZ2K4ZTPN52GK4S7MFXGI4TPNFSA>. Download it today! You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Professor Matthieu J Verstraete Fellow, American Physical Society Chair, Steering Committee, European Theoretical Spectroscopy Facility www.etsf.eu Alumnus Fellow, Young Academy of Europe yacadeuro.org Institute for Theoretical Physics (ITP) Department of Physics Buys Ballot Gebouw/Building, Princetonplein 5, office University of Utrecht, 3584 CC Utrecht ITP Secretariat: +31 30 253 5928 E-mail: ***@***.*** Group web page (Liège): http://www.nanomat.ulg.ac.be/ Nanomat lab, Q-Mat center, Université de Liège Département de Physique, Bat. B5a, 4/50 Allée du 6 août, 19 B-4000 Sart Tilman, Liège Belgium Phone : +32 4 366 90 17 European Theoretical Spectroscopy Facility (ETSF) Mail : ***@***.*** ***@***.*** ***@***.***

[ZhunyunTang]

Dear Matthieu Verstraete,

Thank you for your reply. After conducting further tests, I found that the cause of the error is that the program is unable to properly construct the fourth‑order constraints.

Third-order extraction works flawlessly, even with a very large cutoff (24,606 third‑order parameters, 2,938 third‑order constraints generated successfully). The following is a partial output:
################
total number of equations: 81000
number of FC up 2. order: 2646 --> overdetermination ratio: 30.6
number of FC up 3. order: 24606 --> overdetermination ratio: 3.3
.........
GENERATING CONSTRAINTS
... 62 secondorder constraints (1.30097s)
... 2938 thirdorder constraints (10910.46438s)
... created 3000 constraints (10911.76541s)
... creating fc coefficients 100.0% |========================================| 59.87754s
... least squares solver
... solved second order (98.86606s)
... solved third order (103971.10600s)
... finished solving for forceconstants
..........
###########

For fourth order（In this process, I limited the second- and third-order interactions to the first nearest neighbor）, with far fewer parameters (9,936), the run fails as follows:
###############
total number of equations: 81000
number of FC up 2. order: 270 --> overdetermination ratio: 300.0
number of FC up 3. order: 414 --> overdetermination ratio: 195.7
number of FC up 4. order: 9936 --> overdetermination ratio: 8.2
........
GENERATING CONSTRAINTS
... 62 secondorder constraints (21.98115s)
... 104 thirdorder constraints (14.27390s)
(The code generates 62 second-order and 104 third-order constraints, then immediately segfaults before entering fourth-order constraint generation)
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0 0x14eef84235af in ???
#1 0x55985dfcb90e in ???
#2 0x55985dfe5f7f in ???
#3 0x55985dd1484a in ???
#4 0x55985dbc2ab9 in ???
#5 0x55985dbcff72 in ???
#6 0x14eef840f7e4 in ???
#7 0x55985dbbcb58 in ???
#8 0xffffffffffffffff in ???
###############

Suspected cause:
Given that the problem appears at the constraint-generation step even for a relatively small number of FCs, I suspect a dense-matrix allocation issue. For example, allocate resources based solely on the dimensions of the third-order constraints, leading to out-of-bounds memory access.

Possible way forward:
Is it possible to use a sparse data type (such as the one used in ALAMODE；alamode-code/alamode#213) to avoid dense linear algebra operations?

It should be noted that the above is merely my personal speculation.

Thank you again for your time.

Best regards,
Zhunyun

[mjv500]

Hi Zhunyun Thanks, this is useful sleuthing. Did you manage to recompile with debug options as well? This works give us more to work on, a specific line and subroutine. Did you try to reduce the rc4 progressively to see if it runs with just NN or NNN ? Best Matthieu ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Professor Matthieu J Verstraete Fellow, American Physical Society <https://www.aps.org> Chair, Steering Committee, European Theoretical Spectroscopy Facility <https://www.etsf.eu> Alumnus Fellow, Young Academy of Europe <https://www.yacadeuro.org> <http://yacadeuro.org/> Institute for Theoretical Physics (ITP) Department of Physics Buys Ballot Gebouw/Building, Princetonplein 5, office 7.85 University of Utrecht, 3584 CC Utrecht Phone: +31 30 253 5926 Secretariat: +31 30 253 5929 and ***@***.*** Group web page (Liège): http://www.nanomat.ulg.ac.be/ Nanomat lab, Q-Mat center, Université de Liège Département de Physique, Bat. B5a, office 4/50 Allée du 6 août, 19 B-4000 Sart Tilman, Liège Belgium Phone : +32 4 366 90 17 European Theoretical Spectroscopy Facility <https://www.etsf.eu> (ETSF) Mail : ***@***.*** ***@***.*** ***@***.***

…

On Sat, May 30, 2026, 14:38 Zhunyun Tang ***@***.***> wrote: *ZhunyunTang* left a comment (tdep-developers/tdep#178) <#178 (comment)> Dear Matthieu Verstraete, Thank you for your reply. After conducting further tests, I found that the cause of the error is that the program is unable to properly construct the fourth‑order constraints. Third-order extraction works flawlessly, even with a very large cutoff (24,606 third‑order parameters, 2,938 third‑order constraints generated successfully). The following is a partial output: ################ total number of equations: 81000 number of FC up 2. order: 2646 --> overdetermination ratio: 30.6 number of FC up 3. order: 24606 --> overdetermination ratio: 3.3 ......... GENERATING CONSTRAINTS ... 62 secondorder constraints (1.30097s) ... 2938 thirdorder constraints (10910.46438s) ... created 3000 constraints (10911.76541s) ... creating fc coefficients 100.0% |========================================| 59.87754s ... least squares solver ... solved second order (98.86606s) ... solved third order (103971.10600s) ... finished solving for forceconstants .......... ########### For fourth order（In this process, I limited the second- and third-order interactions to the first nearest neighbor）, with far fewer parameters (9,936), the run fails as follows: ############### total number of equations: 81000 number of FC up 2. order: 270 --> overdetermination ratio: 300.0 number of FC up 3. order: 414 --> overdetermination ratio: 195.7 number of FC up 4. order: 9936 --> overdetermination ratio: 8.2 ........ GENERATING CONSTRAINTS ... 62 secondorder constraints (21.98115s) ... 104 thirdorder constraints (14.27390s) (The code generates 62 second-order and 104 third-order constraints, then immediately segfaults before entering fourth-order constraint generation) Program received signal SIGSEGV: Segmentation fault - invalid memory reference. Backtrace for this error: #0 0x14eef84235af in ??? #1 <#1> 0x55985dfcb90e in ??? #2 <#2> 0x55985dfe5f7f in ??? #3 <#3> 0x55985dd1484a in ??? #4 <#4> 0x55985dbc2ab9 in ??? #5 <#5> 0x55985dbcff72 in ??? #6 <#6> 0x14eef840f7e4 in ??? #7 <#7> 0x55985dbbcb58 in ??? #8 <#8> 0xffffffffffffffff in ??? ############### Suspected cause: Given that the problem appears at the constraint-generation step even for a relatively small number of FCs, I suspect a dense-matrix allocation issue. For example, allocate resources based solely on the dimensions of the third-order constraints, leading to out-of-bounds memory access. Possible way forward: Is it possible to use a sparse data type (such as the one used in ALAMODE； alamode-code/alamode#213 <alamode-code/alamode#213>) to avoid dense linear algebra operations? It should be noted that the above is merely my personal speculation. Thank you again for your time. Best regards, Zhunyun — Reply to this email directly, view it on GitHub <#178?email_source=notifications&email_token=ABWOT32N4H2DDZQYBN7OPS345LI4PA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGI4DGOJSGA42M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4582839209>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABWOT3YTSP6MMTCHXGWMQWT45LI4PAVCNFSM6AAAAACZR4TBI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DKOBSHAZTSMRQHE> . Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS <https://github.com/notifications/mobile/ios/ABWOT32HUJCZK4M7KJQMA7T45LI4PA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGI4DGOJSGA42M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJKTGN5XXIZLSL5UW64Y> and Android <https://github.com/notifications/mobile/android/ABWOT365HGCBBJJLAJSPZWD45LI4PA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGI4DGOJSGA42M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>. Download it today! You are receiving this because you commented.Message ID: ***@***.***>

[ZhunyunTang]

Dear Matthieu Verstraete,

Following your suggestion, I recompiled the code with debug options and ran it again. The error output indeed confirmed our hypothesis – the number of constraints needed for the fourth-order force constants exceeds the preset limit of the array dumeq3. Below is the relevant output：

##########################
At line 1287 of file type_forcemap_constraints.f90
Fortran runtime error: Index '4375' of dimension 2 of array 'dumeq3' outside of expected range (1:4374)

Error termination. Backtrace:
#0 0x556a6ee4f559 in __type_forcemap.type_forcemap_constraints_MOD_lo_fourthorder_asr
at /tdep/src/libolle/type_forcemap_constraints.f90:1287
#1 0x556a6ee70f90 in __type_forcemap_MOD_forceconstant_constraints
at /tdep/tdep-2/1/3/tdep/src/libolle/type_forcemap_constraints.f90:73
#2 0x556a6e91785f in __ifc_solvers_MOD_lo_solve_for_irreducible_ifc
at /tdep/tdep-2/1/3/tdep/src/libolle/ifc_solvers.f90:563
#3 0x556a6e6dbe02 in extract_forceconstants
at /tdep/tdep-2/1/3/tdep/src/extract_forceconstants/main.f90:247
#4 0x556a6e6fe0b2 in main
at /tdep/tdep-2/1/3/tdep/src/extract_forceconstants/main.f90:4

Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.

At line 1287 of file type_forcemap_constraints.f90
Fortran runtime error: Index '4375' of dimension 2 of array 'dumeq3' outside of expected range (1:4374)

Error termination. Backtrace:

mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

Process name: [[9706,1],7]
Exit code: 2

#############################

It appears that when constructing the fourth-order constraints, the required number (likely > 4375) exceeds the current dimension of dumeq3 (which is 4374). Therefore, my question is: can we resolve this by simply increasing the preset size of dumeq3? If so, how would you recommend modifying it? For example, is it enough to change the dimension statement in the source file, or are there other related arrays or logic that need adjustment?

Thank you for your guidance.

Best regards,

Zhunyun

[ZhunyunTang]

Dear Matthieu,

This is good news. By increasing the upper limit of dumeq3, I have successfully resolved the bug. After this modification, the code runs through without any error, and the fourth-order force constants are now extracted successfully.

Best regards,
Zhunyun

[mjv500]

Excellent! I will see about upping the default limit, it may be there based on some benchmarking or just a limit that was previously safe. At any rate there should be an error message explaining the internal limit for the dimensioning. Best Matthieu

…

On Sun, May 31, 2026 at 2:35 PM Zhunyun Tang ***@***.***> wrote: *ZhunyunTang* left a comment (tdep-developers/tdep#178) <#178 (comment)> Dear Matthieu, This is good news. By increasing the upper limit of dumeq3, I have successfully resolved the bug. After this modification, the code runs through without any error, and the fourth-order force constants are now extracted successfully. Best regards, Zhunyun — Reply to this email directly, view it on GitHub <#178?email_source=notifications&email_token=ABWOT3YE5SWX6PMIQ7EMLBT45QRKDA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGY3TAMJWGE22M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4586701615>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABWOT33KBLJQLIYPMD57SD345QRKDAVCNFSM6AAAAACZR4TBI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DKOBWG4YDCNRRGU> . Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS <https://github.com/notifications/mobile/ios/ABWOT324DYG6XI6WF6Q65XD45QRKDA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGY3TAMJWGE22M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJKTGN5XXIZLSL5UW64Y> and Android <https://github.com/notifications/mobile/android/ABWOT32F4P6KLPIKF4GSMVD45QRKDA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTINJYGY3TAMJWGE22M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJLTGN5XXIZLSL5QW4ZDSN5UWI>. Download it today! You are receiving this because you commented.Message ID: ***@***.***>

-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Professor Matthieu J Verstraete Fellow, American Physical Society Chair, Steering Committee, European Theoretical Spectroscopy Facility www.etsf.eu Alumnus Fellow, Young Academy of Europe yacadeuro.org Institute for Theoretical Physics (ITP) Department of Physics Buys Ballot Gebouw/Building, Princetonplein 5, office University of Utrecht, 3584 CC Utrecht ITP Secretariat: +31 30 253 5928 E-mail: ***@***.*** Group web page (Liège): http://www.nanomat.ulg.ac.be/ Nanomat lab, Q-Mat center, Université de Liège Département de Physique, Bat. B5a, 4/50 Allée du 6 août, 19 B-4000 Sart Tilman, Liège Belgium Phone : +32 4 366 90 17 European Theoretical Spectroscopy Facility (ETSF) Mail : ***@***.*** ***@***.*** ***@***.***