search

libAtoms / QUIP

libAtoms/QUIP molecular dynamics framework: https://libatoms.github.io
349 stars 122 forks source link

Phonons at gamma incorrect unless phonon_supercell == phonon_supercell_fine? #36

Closed tdaff closed 8 years ago

tdaff commented 8 years ago

After comparing Daniele Dragoni's commandline for phonons with mine, it seemed like the behaviour phonon_supercell is causing issues. Running the following with recent versions of quippy gives nonzero values for phonons at gamma when phonon_supercell and phonon_supercell_fine differ:

#!/usr/bin/env python

from quippy import phonons
from quippy.potential import Potential
from quippy import bcc1

bulk = bcc1(1.43325007*2, 26)
pot = Potential("IP FS")
pot.minim(bulk, 'cg', 1e-7, 100, do_lat=True)

fine = phonons.Phonon_fine()

steps = 7
sc = 10
sc_fine = 2+(steps*2)

fine.calc(pot, bulk, dx=0.01, phonon_supercell=(sc, sc, sc), 
          phonon_supercell_fine=(sc_fine, sc_fine, sc_fine),
          phonons_path_start=(0.0, 0.0, 0.0), 
          phonons_path_end=(0.5, -0.5, 0.5), phonons_path_steps=steps)

print(sc, fine.q[1], 1000*fine.frequency[1])
print(sc, fine.q[2], 1000*fine.frequency[2])
# (10, FortranArray([ 0.,  0.,  0.]), FortranArray([ 0.68543577,  0.68543577,  0.68543577]))
# (10, FortranArray([ 0.        ,  0.        ,  0.27399202]), FortranArray([ 1.82271397,  1.82271397,  2.40645046]))

sc = sc_fine
fine.calc(pot, bulk, dx=0.01, phonon_supercell=(sc, sc, sc), 
          phonon_supercell_fine=(sc_fine, sc_fine, sc_fine),
          phonons_path_start=(0.0, 0.0, 0.0), 
          phonons_path_end=(0.5, -0.5, 0.5), phonons_path_steps=steps)

print(sc, fine.q[1], 1000*fine.frequency[1])
print(sc, fine.q[2], 1000*fine.frequency[2])
# (16, FortranArray([ 0.,  0.,  0.]), FortranArray([ -2.73119306e-07,  -9.48351971e-08,   2.18828150e-07]))
# (16, FortranArray([ 0.        ,  0.        ,  0.27399202]), FortranArray([ 1.70025326,  1.70025326,  2.37435406]))

It's not a proper bisect, but using a version of quippy compiled in October last year gives the expected results:

(10, FortranArray([ 0.,  0.,  0.]), FortranArray([ -3.96692927e-07,  -2.63923001e-07,   1.37793968e-07]))
(10, FortranArray([ 0.        ,  0.        ,  0.27399202]), FortranArray([ 1.70025326,  1.70025326,  2.37435406]))

(16, FortranArray([ 0.,  0.,  0.]), FortranArray([ -2.65801279e-07,  -1.24993193e-07,   2.24881096e-07]))
(16, FortranArray([ 0.        ,  0.        ,  0.27399202]), FortranArray([ 1.70025326,  1.70025326,  2.37435406]))

Slightly related to this, the documentation for the two parameters is not very clear:

phonon_supercell type=INTEGER dim=3 current_value='1 1 1'
supercell in which to do phonon computation

phonon_supercell_fine type=INTEGER dim=3 current_value='1 1 1'
supercell in which to do phonon computation

Tom

albapa commented 8 years ago

It must be some change I made. I'll take a look.

albapa commented 8 years ago

The problem was the following: we compute forces in phonon_supercell, then we map each atom to atoms in phonon_supercell_fine and copy the forces, as a form of interpolation. I used the criterion of atoms being in the cutoff radius of the potential, but this is obviously wrong as in many cases forces are non-zero until twice the cutoff (or maybe nowhere, in case of Coulombic interactions.) So I changed this so all atoms are mapped and added a warning to catch some silly situations.

I updated the quip command line documentation as well.

schuberm commented 8 years ago

Any ideas why I would see a segmentation fault when I run the above example?

WARNING: Potential_calc: cutoff of Atoms object (-1.00000000000000000) < Potential cutoff (4.40022399999999969) - increasing it now Welcome to minim() space is 12 dimensional Method: Conjugate Gradients Minimization is already converged! Goodbye from minim()

Starting force calculations Displacing atom 1 of 1 Finished force calculations Starting phonon calculations Segmentation fault

I've pulled the latest version.

albapa commented 8 years ago

Hi Samuel,

could you please share your XYZ and potential file (or init args)?

Albert

Sent from my mobile

On 8 Aug 2016 8:11 p.m., "Samuel Huberman" notifications@github.com wrote:

Any ideas why I would see a segmentation fault when I run the above example?

WARNING: Potential_calc: cutoff of Atoms object (-1.00000000000000000) < Potential cutoff (4.40022399999999969) - increasing it now Welcome to minim() space is 12 dimensional Method: Conjugate Gradients Minimization is already converged! Goodbye from minim()

Starting force calculations Displacing atom 1 of 1 Finished force calculations Starting phonon calculations Segmentation fault

I've pulled the latest version.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-238344622, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQiqTCtibbo4rwsRv5duKggbSSZWiks5qd39ggaJpZM4Jb4UU .

schuberm commented 8 years ago

I was simply running the example @tdaff gave in the first comment (python example.py), which I believe references the ip.parms.FS.xml in ${QUIP_ROOT}/share/Parameters. Have I misunderstood how this example should be used?

jameskermode commented 8 years ago

If you are using the intel compiler suite it could be a stack size issue when reading the XML file: try ulimit -S unlimited.

albapa commented 8 years ago

OK, I understand now. I will take a look a bit later.

After updating, did you do "make clean" and then "make install-quippy"?

Albert

Sent from my mobile

On 8 Aug 2016 8:56 p.m., "Samuel Huberman" notifications@github.com wrote:

I was simply running the example @tdaff https://github.com/tdaff gave in the first comment (python example.py), which I believe references the ip.parms.FS.xml in ${QUIP_ROOT}/share/Parameters. Have I misunderstood how this example should be used?

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-238357080, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQqk4FzFY6cqkcDQQr1syF6MqJ8pHks5qd4n5gaJpZM4Jb4UU .

schuberm commented 8 years ago

Ran make deepclean, make, make quippy with export QUIP_ARCH=linux_x86_64_gfortran, still see the seg. fault.

albapa commented 8 years ago

Weird - I tried on both linux and mac and the script runs fine.

What's in your PYTHONPATH? In other words, is your python environment using the latest quippy that you have just compiled?

Albert

On 8 August 2016 at 21:35, Samuel Huberman notifications@github.com wrote:

Ran make deepclean, make, make quippy with export QUIP_ARCH=linux_x8664 gfortran, still see the seg. fault.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-238368521, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQlU1Y8uVu8sXHBxr73EfobkPWfyHks5qd5MugaJpZM4Jb4UU .

schuberm commented 8 years ago

Ya it seems weird.

Here's what my .bashrc looks like:

export PYTHONPATH=/home/schuberm/installs/numpy-1.11.1rc1/lib/python2.7/site-packages:$PYTHONPATH
export PYTHONPATH=/home/schuberm/installs/scipy-0.18.0rc1//lib/python2.7/site-packages:$PYTHONPATH
export PYTHONPATH=/home/schuberm/installs/QUIP/quippy/build/linux_x86_64_gfortran/lib.linux-x86_64-2.7:$PYTHONPATH

Here is my Makefile.inc

# Place to override setting elsewhere, in particular things set in Makefile.linux_x86_64_gfortran
# look in /home/schuberm/installs/QUIP/arch/Makefile.linux_x86_64_gfortran for defaults set by arch
# OPTIM=
# COPTIM=
# DEBUG=
# DEBUG=-O0 -g -DDUMP_CORE_ON_ABORT -DDEBUG -fbounds-check
# CDEBUG=
MATH_LINKOPTS=-L/cm/shared/apps/lapack/gcc/64/3.4.2 -llapack -L/cm/shared/apps/openblas/dynamic/0.2.6/lib -lopenblas
HAVE_NETCDF=0
USE_MAKEDEP=1
MAKEDEP=makedep.py
HAVE_MDCORE=0
HAVE_ASAP=0
HAVE_KIM=0
HAVE_CP2K=0
HAVE_VASP=1
HAVE_TB=1
HAVE_PRECON=1
HAVE_FX=0
HAVE_GAP=0
HAVE_GAP_FILLER=0
HAVE_QR=1
HAVE_LOTF=0
HAVE_ONIOM=0
HAVE_LOCAL_E_MIX=0
HAVE_QC=0
EXTRA_LINKOPTS=
HAVE_CGAL=0
SIZEOF_FORTRAN_T= 2
HAVE_THIRDPARTY=0
HAVE_METIS=0
HAVE_SCME=0
HAVE_MTP=0
HAVE_LMTO_TBE=0
albapa commented 8 years ago

Which version of gfortran are you using?

On 8 August 2016 at 23:41, Samuel Huberman notifications@github.com wrote:

Ya it's seems weird.

Here's what my .bashrc looks like:

export PYTHONPATH=/home/schuberm/installs/numpy-1.11.1rc1/lib/python2.7/site-packages:$PYTHONPATH export PYTHONPATH=/home/schuberm/installs/scipy-0.18.0rc1//lib/python2.7/site-packages:$PYTHONPATH export PYTHONPATH=/home/schuberm/installs/QUIP/quippy/build/linux_x86_64_gfortran/lib.linux-x86_64-2.7:$PYTHONPATH

Here is my Makefile.inc

Place to override setting elsewhere, in particular things set in Makefile.linux_x86_64_gfortran

look in /home/schuberm/installs/QUIP/arch/Makefile.linux_x86_64_gfortran for defaults set by arch

OPTIM=

COPTIM=

DEBUG=

DEBUG=-O0 -g -DDUMP_CORE_ON_ABORT -DDEBUG -fbounds-check

CDEBUG=

MATH_LINKOPTS=-L/cm/shared/apps/lapack/gcc/64/3.4.2 -llapack -L/cm/shared/apps/openblas/dynamic/0.2.6/lib -lopenblas HAVE_NETCDF=0 USE_MAKEDEP=1 MAKEDEP=makedep.py HAVE_MDCORE=0 HAVE_ASAP=0 HAVE_KIM=0 HAVE_CP2K=0 HAVE_VASP=1 HAVE_TB=1 HAVE_PRECON=1 HAVE_FX=0 HAVE_GAP=0 HAVE_GAP_FILLER=0 HAVE_QR=1 HAVE_LOTF=0 HAVE_ONIOM=0 HAVE_LOCAL_E_MIX=0 HAVE_QC=0 EXTRA_LINKOPTS= HAVE_CGAL=0 SIZEOF_FORTRAN_T= 2 HAVE_THIRDPARTY=0 HAVE_METIS=0 HAVE_SCME=0 HAVE_MTP=0 HAVE_LMTO_TBE=0

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-238399992, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQg-YLP6LbRvpQLr6AaAFeatd-Yz7ks5qd7B_gaJpZM4Jb4UU .

schuberm commented 8 years ago 
gfortran -v
Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/cm/shared/apps/gcc/4.8.1/libexec/gcc/x86_64-unknown-linux-gnu/4.8.1/lto-wrapper
Target: x86_64-unknown-linux-gnu
Configured with: ../gcc-4.8.1/configure --prefix=/cm/shared/apps/gcc/4.8.1 --enable-languages=c,c++,fortran --with-gmp-include=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-gmp --with-gmp-lib=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-gmp --with-mpc-include=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpc/src --with-mpc-lib=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpc/src/.libs --with-mpfr-include=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpfr/src --with-mpfr-lib=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpfr/src/.libs
Thread model: posix
gcc version 4.8.1 (GCC)
albapa commented 8 years ago

Hmm, I tried with both 4.8.0 and 4.8.5 versions of gfortran (those are the ones I could access easily) and Tom's script finishes without an error in both cases.

I also ran it with debug options enabled, without any problems.

Could you perhaps try with another compiler and/or compiling with debug options? I am sorry for not being too helpful...

Albert

On 9 August 2016 at 00:09, Samuel Huberman notifications@github.com wrote:

gfortran -v Using built-in specs. COLLECT_GCC=gfortran COLLECT_LTO_WRAPPER=/cm/shared/apps/gcc/4.8.1/libexec/gcc/x86_64-unknown-linux-gnu/4.8.1/lto-wrapper Target: x86_64-unknown-linux-gnu Configured with: ../gcc-4.8.1/configure --prefix=/cm/shared/apps/gcc/4.8.1 --enable-languages=c,c++,fortran --with-gmp-include=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-gmp --with-gmp-lib=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-gmp --with-mpc-include=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpc/src --with-mpc-lib=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpc/src/.libs --with-mpfr-include=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpfr/src --with-mpfr-lib=/usr/src/redhat/BUILD/gcc-4.8.1-obj/../gcc-4.8.1/our-mpfr/src/.libs Thread model: posix gcc version 4.8.1 (GCC)

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-238405256, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQtAla3dLi1HsT8ABrdqR_W7zdrKgks5qd7czgaJpZM4Jb4UU .

jameskermode commented 8 years ago

Hi Samuel, I agree with Albert, the next step is re-compiling with the DEBUG="-O0 -g" ... and OPTIM= lines in the Makefile.inc uncommented. If you then run the script from within the gdb debugger (gdb python, then run example.py) you should get a stack trace when the seg fault happens, which will help us figure out where it's crashing.

jameskermode commented 8 years ago

(I just checked, the script runs fine with gfortran 4.9 on my Mac)

schuberm commented 8 years ago

Looks like openblas is the culprit, output from gdb

Welcome to minim()
space is 12 dimensional
Method: Conjugate Gradients
Minimization is already converged!
Goodbye from minim()

Starting force calculations
Displacing atom 1 of 1
Finished force calculations
Starting phonon calculations

Program received signal SIGSEGV, Segmentation fault.
0x00002aaab839c006 in ?? () from /cm/shared/apps/openblas/dynamic/0.2.6/lib/libopenblas.so.0
Missing separate debuginfos, use: debuginfo-install glibc-2.12-1.107.el6.x86_64 keyutils-libs-1.4-4.el6.x86_64 krb5-libs-1.10.3-10.el6.x86_64 libcom_err-1.41.12-14.el6.x86_64 libselinux-2.0.94-5.3.el6.x86_64 openssl-1.0.0-27.el6.x86_64 zlib-1.2.3-29.el6.x86_64
jameskermode commented 8 years ago

Can you get a full stack trace? (bt command in gdb)

schuberm commented 8 years ago 
#0  0x00002aaab839bf2d in ?? () from /cm/shared/apps/openblas/dynamic/0.2.6/lib/libopenblas.so.0
#1  0x00002aaab8063d21 in zdotc_ () from /cm/shared/apps/openblas/dynamic/0.2.6/lib/libopenblas.so.0
#2  0x00002aaab7b3d244 in zhetd2_ () from /cm/shared/apps/lapack/gcc/64/3.4.2/liblapack.so
#3  0x00002aaab7b3fb8e in zhetrd_ () from /cm/shared/apps/lapack/gcc/64/3.4.2/liblapack.so
#4  0x00002aaab7b34836 in zheev_ () from /cm/shared/apps/lapack/gcc/64/3.4.2/liblapack.so
#5  0x00002aaab61af885 in __linearalgebra_module_MOD_matrix_z_diagonalise ()
   from /home/schuberm/installs/QUIP/quippy/build/linux_x86_64_gfortran/lib.linux-x86_64-2.7/quippy/_quippy.so
#6  0x00002aaab5d97bf2 in __phonons_module_MOD_phonon_fine_calc ()
   from /home/schuberm/installs/QUIP/quippy/build/linux_x86_64_gfortran/lib.linux-x86_64-2.7/quippy/_quippy.so
#7  0x00002aaab5c65417 in f2py_rout__quippy_qp_phonon_fine_calc ()
    at build/linux_x86_64_gfortran/src.linux-x86_64-2.7/quippy/_quippymodule.c:168454
#8  0x0000000000422e1a in PyObject_Call ()
#9  0x00000000004b26d2 in PyEval_EvalFrameEx () at Python/ceval.c:4663
#10 0x00000000004ba108 in PyEval_EvalCodeEx () at Python/ceval.c:3582
#11 0x000000000052ef9f in function_call ()
#12 0x0000000000422e1a in PyObject_Call ()
#13 0x00000000004b26d2 in PyEval_EvalFrameEx () at Python/ceval.c:4663
#14 0x00000000004ba108 in PyEval_EvalCodeEx () at Python/ceval.c:3582
#15 0x00000000004b8c24 in PyEval_EvalFrameEx () at Python/ceval.c:4446
#16 0x00000000004ba108 in PyEval_EvalCodeEx () at Python/ceval.c:3582
#17 0x00000000004ba232 in PyEval_EvalCode () at Python/ceval.c:669
#18 0x00000000004e77dd in PyRun_FileExFlags () at Python/pythonrun.c:1370
#19 0x00000000004e9061 in PyRun_SimpleFileExFlags () at Python/pythonrun.c:948
#20 0x0000000000415a7d in Py_Main ()
#21 0x0000003f73c1ecdd in __libc_start_main () from /lib64/libc.so.6
#22 0x0000000000414bf5 in _start ()
jameskermode commented 8 years ago

The matrix diagonalisation is failing - most likely a problem with your OpenBLAS installation. Can you try with a different version of BLAS and/or LAPACK?

schuberm commented 8 years ago

Recompiled with libblas.so, runs without seg. fault.

tdaff commented 8 years ago

@schuberm I've managed to produce some segfaults with OpenBLAS that I tracked down to threading issues. If you wanted to try it again (it's faster than reference BLAS), you could try export OPENBLAS_NUM_THREADS=1 to stop it threading or recompile OpenBLAS with threading off or OpenMP on.

schuberm commented 8 years ago

Sorry to keep posting on this issue, but I have some more questions about the phonon calculation. Where would be an appropriate place to ask?

albapa commented 8 years ago

I think asking here is fine but you can also write to me directly.

Albert

Sent from my mobile

On 25 Aug 2016 6:07 p.m., "Samuel Huberman" notifications@github.com wrote:

Sorry to keep posting on this issue, but I have some more questions about the phonon calculation. Where would be an appropriate place to ask?

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-242466649, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQp09dlsbmTNuoFkgwi-4jH5YHnbrks5qjcvlgaJpZM4Jb4UU .

schuberm commented 8 years ago

I'm trying to understand the output of a Phonon_fine.calc call. So in the above script

fine = phonons.Phonon_fine()

steps = 7
sc = 10
sc_fine = 2+(steps*2)

fine.calc(pot, bulk, dx=0.01, phonon_supercell=(sc, sc, sc), 
          phonon_supercell_fine=(sc_fine, sc_fine, sc_fine),
          phonons_path_start=(0.0, 0.0, 0.0), 
          phonons_path_end=(0.5, -0.5, 0.5), phonons_path_steps=steps)

The descriptor fine.eigenvector contains the eigenvectors in a [3, number of atoms, 3*number of atoms, steps+2] array. So a given eigenvector of a single mode is in the [:,:, x,y] slice of the array?

albapa commented 8 years ago

Yes, that's correct, see lines 101-103 in src/Utils/phonons.f95:

  allocate(this%q(3,this%n_qvectors))
  allocate(this%frequency(this%n_modes,this%n_qvectors))
  allocate(this%eigenvector(3,n_modes/3,n_modes,n_qvectors))

and n_modes = 3 * number of atoms.

I hope this helps.

Albert

On 25 August 2016 at 19:20, Samuel Huberman notifications@github.com wrote:

I'm trying to understand the output of a Phonon_fine.calc call. So in the above script

fine = phonons.Phonon_fine()

steps = 7 sc = 10 sc_fine = 2+(steps*2)

fine.calc(pot, bulk, dx=0.01, phonon_supercell=(sc, sc, sc), phonon_supercell_fine=(sc_fine, sc_fine, sc_fine), phonons_path_start=(0.0, 0.0, 0.0), phonons_path_end=(0.5, -0.5, 0.5), phonons_path_steps=steps)

The descriptor fine.eigenvector contains the eigenvectors in a [3, number of atoms, 3*number of atoms, steps+2] array. So a given eigenvector of a single mode is in the [:,:, x,y] slice of the array?

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/libAtoms/QUIP/issues/36#issuecomment-242490195, or mute the thread https://github.com/notifications/unsubscribe-auth/AEKWQmicDdjTVLdCSWyJbsLeq8TXkLcWks5qjd0AgaJpZM4Jb4UU .

schuberm commented 8 years ago

Thanks. I'm also wondering about line 370 src/Utils/phonons.f95:

this%eigenvector(:,:,:,k) = real( reshape( evecs, (/3, at_in%N, 3*at_in%N /) ), dp )

Why do we take the real component here? I believe the normal modes can be complex, and only at the gamma point do the components take purely real values.

albapa commented 8 years ago

Yes, you're right, this was carried from the past and never got updated. I changed this so if you do an update the new version should return the proper modes. Thanks for pointing this out.