diff --git a/lib/README b/lib/README index 192589c8b0..ffb3e3f446 100644 --- a/lib/README +++ b/lib/README @@ -11,7 +11,9 @@ The libraries included in the LAMMPS distribution are the following: atc atomistic-to-continuum methods from Reese Jones, Jeremy Templeton, Jon Zimmerman (Sandia) -gpu graphical processor (GPU) routines, currently NVIDIA specific, +awpmd antisymmetrized wave packet molecular dynamics + from Ilya Valuev (JIHT RAS) +gpu graphical processor (GPU) routines, currently NVIDIA specific from Mike Brown (Sandia) poems POEMS rigid-body integration package from RPI diff --git a/lib/awpmd/Makefile.openmpi b/lib/awpmd/Makefile.openmpi new file mode 100644 index 0000000000..d70754ef16 --- /dev/null +++ b/lib/awpmd/Makefile.openmpi @@ -0,0 +1,64 @@ +SHELL = /bin/sh + +# ------ FILES ------ + +SRC = \ + ivutils/src/logexc.cpp \ + systems/interact/TCP/wpmd.cpp \ + systems/interact/TCP/wpmd_split.cpp + +INC = \ + cerf.h \ + cerf2.h \ + cerf_octave.h \ + cvector_3.h \ + lapack_inter.h \ + logexc.h \ + pairhash.h \ + refobj.h \ + tcpdefs.h \ + vector_3.h \ + wavepacket.h \ + wpmd.h \ + wpmd_split.h + +# ------ DEFINITIONS ------ + +LIB = libawpmd.a +OBJ = $(SRC:.cpp=.o) + +# ------ SETTINGS ------ + +# include any MPI settings needed for the ATC library to build with +# the same MPI library that LAMMPS is built with + +CC = mpic++ +CCFLAGS = -O -Isystems/interact/TCP/ -Isystems/interact -Iivutils/include +ARCHIVE = ar +ARCHFLAG = -rc +DEPFLAGS = -M +#LINK = +#LINKFLAGS = +USRLIB = +SYSLIB = + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + +# ------ COMPILE RULES ------ + +%.o:%.cpp + $(CC) $(CCFLAGS) -c $< -o $@ +%.d:%.cpp + $(CC) $(CCFLAGS) $(DEPFLAGS) $< > $@ + +# ------ DEPENDENCIES ------ + +DEPENDS = $(OBJ:.o=.d) + +# ------ CLEAN ------ + +clean: + rm *.d *~ $(OBJ) $(LIB) diff --git a/lib/awpmd/README b/lib/awpmd/README new file mode 100644 index 0000000000..0806a31b39 --- /dev/null +++ b/lib/awpmd/README @@ -0,0 +1,29 @@ +AWPMD (Antisymmetrized Wave Packet Molecular Dynamics) library + +Ilya Valuev, Igor Morozov, JIHT RAS +valuev at physik.hu-berlin.de +June 2011 + +-------------- + +This is version 0.9 of the AWPMD library taken from JIHT GridMD project. +It contains interface to calculate electronic and electron-ion Hamiltonian, +norm matrix and forces for AWPMD method. + + +This library must be built with a C++ compiler, before LAMMPS is +built, so LAMMPS can link against it. + +Build the library using one of the provided Makefiles or creating your +own, specific to your compiler and system. For example: + +make -f Makefile.openmpi++ + +If the build is successful, you should end up with a libawpmd.a file. + +-------------- + +AWPMD is an open source program distributed under the terms +of wxWidgets Library License (see license directory for details). + + diff --git a/lib/awpmd/ivutils/include/cerf.h b/lib/awpmd/ivutils/include/cerf.h new file mode 100644 index 0000000000..207fe8f459 --- /dev/null +++ b/lib/awpmd/ivutils/include/cerf.h @@ -0,0 +1,245 @@ +# ifndef CERF_H +# define CERF_H + +# include +# include +/* + +Copyright (C) 1998, 1999 John Smith + +This file is part of Octave. +Or it might be one day.... + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Put together by John Smith john at arrows dot demon dot co dot uk, +// using ideas by others. +// +// Calculate erf(z) for complex z. +// Three methods are implemented; which one is used depends on z. +// +// The code includes some hard coded constants that are intended to +// give about 14 decimal places of accuracy. This is appropriate for +// 64-bit floating point numbers. +// +// Oct 1999: Fixed a typo that in +// const Complex cerf_continued_fraction( const Complex z ) +// that caused erroneous answers for erf(z) where real(z) negative +// + + +//#include +//#include + +//#include "f77-fcn.h" + +// +// Abramowitz and Stegun: (eqn: 7.1.14) gives this continued +// fraction for erfc(z) +// +// erfc(z) = sqrt(pi).exp(-z^2). 1 1/2 1 3/2 2 5/2 +// --- --- --- --- --- --- ... +// z + z + z + z + z + z + +// +// This is evaluated using Lentz's method, as described in the narative +// of Numerical Recipes in C. +// +// The continued fraction is true providing real(z)>0. In practice we +// like real(z) to be significantly greater than 0, say greater than 0.5. +// +template< class Complex> +const Complex cerfc_continued_fraction( const Complex z ) +{ + double tiny = 1e-20 ; // a small number, large enough to calculate 1/tiny + double eps = 1e-15 ; // large enough so that 1.0+eps > 1.0, when using + // the floating point arithmetic + // + // first calculate z+ 1/2 1 + // --- --- ... + // z + z + + Complex f(z) ; + Complex C(f) ; + Complex D(0.0) ; + Complex delta ; + double a ; + + a = 0.0 ; + do + { + a = a + 0.5 ; + D = z + a*D ; + C = z + a/C ; + + if (D.real() == 0.0 && D.imag() == 0.0) + D = tiny ; + + D = 1.0 / D ; + + delta = (C * D) ; + + f = f * delta ; + + } while (abs(1.0-delta) > eps ) ; + + // + // Do the first term of the continued fraction + // + f = 1.0 / f ; + + // + // and do the final scaling + // + f = f * exp(-z*z)/ sqrt(M_PI) ; + + return f ; +} + +template< class Complex> +const Complex cerf_continued_fraction( const Complex z ) +{ + // warning("cerf_continued_fraction:"); + if (z.real() > 0) + return 1.0 - cerfc_continued_fraction( z ) ; + else + return -1.0 + cerfc_continued_fraction( -z ) ; +} + +// +// Abramawitz and Stegun, Eqn. 7.1.5 gives a series for erf(z) +// good for all z, but converges faster for smallish abs(z), say abs(z)<2. +// +template< class Complex> +const Complex cerf_series( const Complex z ) +{ + double tiny = 1e-20 ; // a small number compared with 1. + // warning("cerf_series:"); + Complex sum(0.0) ; + Complex term(z) ; + Complex z2(z*z) ; + + for (int n=0; n<3 || abs(term) > abs(sum)*tiny; n++) + { + sum = sum + term / (2*n+1) ; + term = -term * z2 / (n+1) ; + } + + return sum * 2.0 / sqrt(M_PI) ; +} + +// +// Numerical Recipes quotes a formula due to Rybicki for evaluating +// Dawson's Integral: +// +// exp(-x^2) integral exp(t^2).dt = 1/sqrt(pi) lim sum exp(-(z-n.h)^2) / n +// 0 to x h->0 n odd +// +// This can be adapted to erf(z). +// +template< class Complex> +const Complex cerf_rybicki( const Complex z ) +{ + // warning("cerf_rybicki:"); + double h = 0.2 ; // numerical experiment suggests this is small enough + + // + // choose an even n0, and then shift z->z-n0.h and n->n-h. + // n0 is chosen so that real((z-n0.h)^2) is as small as possible. + // + int n0 = 2*(int) (floor( z.imag()/(2*h) + 0.5 )) ; + + Complex z0( 0.0, n0*h ) ; + Complex zp(z-z0) ; + Complex sum(0.0,0.0) ; + // + // limits of sum chosen so that the end sums of the sum are + // fairly small. In this case exp(-(35.h)^2)=5e-22 + // + // + for (int np=-35; np<=35; np+=2) + { + Complex t( zp.real(), zp.imag()-np*h) ; + Complex b( exp(t*t) / (np+n0) ) ; + sum += b ; + } + + sum = sum * 2 * exp(-z*z) / M_PI ; + + return Complex(-sum.imag(), sum.real()) ; +} + +template< class Complex> +const Complex cerf( const Complex z ) +{ + // + // Use the method appropriate to size of z - + // there probably ought to be an extra option for NaN z, or infinite z + // + // + if (abs(z) < 2.0) + return cerf_series( z ) ; + else if (abs(z.real()) < 0.5) + return cerf_rybicki( z ) ; + else + return cerf_continued_fraction( z ) ; +} + +// +// Footnote: +// +// Using the definitions from Abramowitz and Stegun (7.3.1, 7.3.2) +// The fresnel intgerals defined as: +// +// / t=x +// C(x) = | cos(pi/2 t^2) dt +// / +// t=0 +// +// and +// / t=x +// S(x) = | sin(pi/2 t^2) dt +// / +// t=0 +// +// These can be derived from erf(x) using 7.3.22 +// +// C(z) +iS(z) = (1+i) erf( sqrt(pi)/2 (1-i) z ) +// ----- +// 2 +// +// -------------------------------------------------------------------------- +// Some test examples - +// comparative data taken from Abramowitz and Stegun table 7.9. +// Table 7.9 tabulates w(z), where w(z) = exp(-z*z) erfc(iz) +// I have copied twelve values of w(z) from the table, and separately +// calculated them using this code. The results are identical. +// +// x y Abramowitz & Stegun | Octave Calculations +// w(x+iy) | w(x+iy) cerf ( i.(x+iy)) +// 0.2 0.2 0.783538+0.157403i | 0.783538 +0.157403 0.23154672 -0.219516 +// 0.2 0.7 0.515991+0.077275i | 0.515991 +0.077275 0.69741968 -0.138277 +// 0.2 1.7 0.289309+0.027154i | 0.289309 +0.027154 0.98797507 -0.011744 +// 0.2 2.7 0.196050+0.013002i | 0.196050 +0.013002 0.99994252 -0.000127 +// 1.2 0.2 0.270928+0.469488i | 0.270928 +0.469488 0.90465623 -2.196064 +// 1.2 0.7 0.280740+0.291851i | 0.280740 +0.291851 1.82926135 -0.639343 +// 1.2 1.7 0.222436+0.129684i | 0.222436 +0.129684 1.00630308 +0.060067 +// 1.2 2.7 0.170538+0.068617i | 0.170538 +0.068617 0.99955699 -0.000290 +// 2.2 0.2 0.041927+0.287771i | 0.041927 +0.287771 24.70460755-26.205981 +// 2.2 0.7 0.099943+0.242947i | 0.099943 +0.242947 9.88734713+18.310797 +// 2.2 1.7 0.135021+0.153161i | 0.135021 +0.153161 1.65541359 -1.276707 +// 2.2 2.7 0.127900+0.096330i | 0.127900 +0.096330 0.98619434 +0.000564 + +# endif \ No newline at end of file diff --git a/lib/awpmd/ivutils/include/cvector_3.h b/lib/awpmd/ivutils/include/cvector_3.h new file mode 100644 index 0000000000..be2296275d --- /dev/null +++ b/lib/awpmd/ivutils/include/cvector_3.h @@ -0,0 +1,75 @@ +/*s*************************************************************************** + * + * Copyright (c), Ilya Valuev 2005 All Rights Reserved. + * + * Author : Ilya Valuev, MIPT, Moscow, Russia + * + * Project : GridMD, ivutils + * + * + * + *****************************************************************************/ +/*r @file vector_3.h @brief работа с трехмерными комплексными векторами +*/ + +# ifndef CVECTOR_3A_H +# define CVECTOR_3A_H + +# include +# include +# include "vector_3.h" + +using namespace std; + +typedef complex cdouble; +typedef Vector_Nt cVector_3; + +//------------------------------------------------------ +// Overloads for cdouble +//------------------------------------------------------ + +inline cdouble operator*(int a, const cdouble &b){ + return ((double)a)*b; +} +inline cdouble operator*(const cdouble &b,int a){ + return a*b; +} +inline cdouble operator/(const cdouble &b,int a){ + return (1./a)*b; +} + +inline cdouble operator/(int a, const cdouble &b){ + return (a)*(1./b); +} + +//------------------------------------------------------ +// Overloads for cVector_3 +//------------------------------------------------------ + +inline cVector_3 operator*(const cdouble &a, const Vector_3 &v){ + return cVector_3(a*v[0], a*v[1], a*v[2]); // a*cVector_3(v); +} + +inline cVector_3 operator*(const Vector_3 &v, const cdouble &a){ + return a*v; +} + +inline Vector_3 real(const cVector_3 &cv){ + return Vector_3(cv[0].real(), cv[1].real(), cv[2].real()); +} + +inline Vector_3 imag(const cVector_3 &cv){ + return Vector_3(cv[0].imag(), cv[1].imag(), cv[2].imag()); +} + +inline cVector_3 conj(const cVector_3 &cv){ + return cVector_3(conj(cv[0]), conj(cv[1]), conj(cv[2])); +} + +inline cVector_3 rcell1(cVector_3& cv, Vector_3 &cell, int flags=0xffff) { + return cVector_3( real(cv).rcell1(cell, flags) ) + cdouble(0,1)*imag(cv); +} + + +# endif // __CVECTOR_3A_H + diff --git a/lib/awpmd/ivutils/include/erf.h b/lib/awpmd/ivutils/include/erf.h new file mode 100644 index 0000000000..b1cdbcb9dc --- /dev/null +++ b/lib/awpmd/ivutils/include/erf.h @@ -0,0 +1,19 @@ +# ifndef ERF_H +# define ERF_H + +# ifdef _WIN32 + +# ifdef __cplusplus +extern "C" { +# endif + +double erf(double x); +double erfc(double x); + +# ifdef __cplusplus +} +# endif + +# endif + +# endif diff --git a/lib/awpmd/ivutils/include/lapack_inter.h b/lib/awpmd/ivutils/include/lapack_inter.h new file mode 100644 index 0000000000..743c5e6059 --- /dev/null +++ b/lib/awpmd/ivutils/include/lapack_inter.h @@ -0,0 +1,53 @@ +// Interface for LAPACK function + +# ifndef LAPACK_INTER_H +# define LAPACK_INTER_H + +#include +typedef int lapack_int; +typedef complex lapack_complex_float; +typedef complex lapack_complex_double; + +#ifdef _WIN32 + + //#define MKL_Complex8 lapack_complex_float + //#define MKL_Complex16 lapack_complex_double + #include "mkl.h" + + inline void ZPPTRF( char* uplo, const lapack_int* n, lapack_complex_double* ap, lapack_int* info ) { + ZPPTRF(uplo, (int*)n, (MKL_Complex16*)ap, (int*)info); + } + inline void ZPPTRI( char* uplo, const lapack_int* n, lapack_complex_double* ap, lapack_int* info ){ + ZPPTRI(uplo, (int*)n, (MKL_Complex16*)ap, (int*)info); + } + +#else + + #define DGETRF dgetrf_ + #define DGETRS dgetrs_ + #define DGETRI dgetri_ + #define ZPPTRF zpptrf_ + #define ZPPTRI zpptri_ + + #ifdef __cplusplus + extern "C" { + #endif /* __cplusplus */ + void dgetrf_( const lapack_int* m, const lapack_int* n, double* a, const lapack_int* lda, + lapack_int* ipiv, lapack_int* info ); + void dgetrs_( const char* trans, const lapack_int* n, const lapack_int* nrhs, + const double* a, const lapack_int* lda, const lapack_int* ipiv, + double* b, const lapack_int* ldb, lapack_int* info ); + void dgetri_( const lapack_int* n, double* a, const lapack_int* lda, + const lapack_int* ipiv, double* work, const lapack_int* lwork, + lapack_int* info ); + void zpptrf_( const char* uplo, const lapack_int* n, lapack_complex_double* ap, + lapack_int* info ); + void zpptri_( const char* uplo, const lapack_int* n, lapack_complex_double* ap, + lapack_int* info ); + #ifdef __cplusplus + } + #endif /* __cplusplus */ + +#endif + +#endif /* lapack_intER_H */ diff --git a/lib/awpmd/ivutils/include/logexc.h b/lib/awpmd/ivutils/include/logexc.h new file mode 100644 index 0000000000..8d6bc50fe2 --- /dev/null +++ b/lib/awpmd/ivutils/include/logexc.h @@ -0,0 +1,315 @@ +# ifndef LOGEXC_H +# define LOGEXC_H + +#include +#include +#include +#include + +# ifndef _WIN32 +# include +# endif + + +/// this specifies whether to put file/line info in error messages +# ifndef LINFO +# define LINFO 1 +# endif + +using namespace std; + +/// verbosity levels +enum vbLEVELS{ + vblNONE = 0, ///< completely silent + vblFATAL = 0x1, ///< report fatal errors + vblOERR = 0x2, ///< report other errors + vblERR = vblFATAL|vblOERR, ///< report all errors + vblWARN = 0x4, ///< report warnings + vblALLBAD = vblWARN|vblERR, ///< report all errors and warnings + vblMESS1 = 0x8, ///< report messages of level 1 (important) + vblMESS2 = 0x10, ///< report messages of level 2 (less important) + vblMESS3 = 0x20, ///< report messages of level 3 (not important) + vblMESS4 = 0x40, ///< report messages of level 4 (anything possible) + vblALLMESS = vblMESS1|vblMESS2|vblMESS3|vblMESS4, ///< report all messages + vblPROGR = 0x80, ///< report progress + vblALL = 0xffff +}; + + + +/// traits structure to deduce exception level and name from exception class +/// by default all exceptions have vblFATAL level +template +struct log_exception_traits{ + /// exeption level according to the vbLEVELS + static int level(const exc_t &signal){ return vblFATAL; } + /// the string name of exception category + static string name(const exc_t &signal){ return typeid(exc_t).name();} + /// adds some more explanations to the description + /// default behaviour: nothing done + static exc_t add_words(const exc_t &orig, const char *words){ + return orig; + } +}; + + + + +/// integer exceptions have the level equal to their value +template<> +struct log_exception_traits{ + /// exeption level according to the vbLEVELS + static int level(const int &signal){ return signal; } + /// the string name of exception category + static string name(const int &signal){ + if(signal&vblFATAL) + return "fatal error"; + else if(signal&vblOERR) + return "error"; + else if(signal&vblWARN) + return "warning"; + else + return ""; + /* + else if(signal&vblALLMESS) + return "message"; + else if(signal&vblPROGR) + return "progress report"; + else + return "integer exception";*/ + } + /// default behaviour: nothing done + static int add_words(const int &orig, const char *words){ + return orig; + } +}; + +/// vbLEVELS exceptions act as integers +template<> +struct log_exception_traits{ + static int level(const enum vbLEVELS &signal){ return log_exception_traits::level(signal); } + static string name(const enum vbLEVELS &signal){ return log_exception_traits::name(signal); } + static enum vbLEVELS add_words(const enum vbLEVELS &orig, const char *words){ + return orig; + } +}; + + + + +/// Logger class to control (computational) function behaviour when something requiring user attention has happened. +/// message(signal,errcode, text) is used to either throw an exception or return errorcode +/// At first, the the level of error is determined via log_exception_traits<>::level(signal) +/// For integer (enum) signals the level is the signal itself. +/// Then text is printed, if signal level is listed in output levels or (or in extra outlevels, if they are set) +/// via log_text() function. +/// If level has vblERR bit, the behaviour is controlled by the flag specified in set_throw(flag): +/// flag=0: nothing done; +/// flag=1: calls add_words() for signal and throws signal; +/// flag=2: throws pair<>(errcode, text); +/// flag=3: throws errcode. +/// Then, if the level is listed in stop_levels (or in extra stop levels, if they are set), the program is aborted, +/// otherwise errcode is returned; +/// The function set_levels(out_levels,stop_levels) is used to specify bitflags for the levels which +/// require message output or/and program termination. Stop level has effect only when exceptions are not thrown. +/// The function extra_levels(eout_levels,estop_levels) is used to temporarily set the corresponding levels, +/// they are unset (the original levels are restored) by calling extra_levels(0,0). +class message_logger { + // global message is a friend + template + friend int message(const exc_t &signal, int errcode, const char *what, ...); +protected: + string descriptor; + int throw_ex; + int outlevel, eoutlevel; + int stoplevel, estoplevel; + + static message_logger *glogp; + /// used to restore the previous global logger + message_logger *prev, *next; +public: + + message_logger(const string &descriptor_="", int out_level=vblALLBAD|vblMESS1, + int stop_level=vblFATAL, int throw_exceptions=0, int use_globally=0) + :descriptor(descriptor_),prev(NULL),next(NULL){ + set_throw(throw_exceptions); + set_levels(out_level,stop_level); + extra_levels(0,0); + if(use_globally){ + set_global(true); + } + } + + /// returns a reference to global logger + /// if not set, links with default message_logger + static message_logger &global(); + + /// sets/unsets this logger as the global logger + int set_global(bool set){ + if(set){ + if(prev) // already set + return -1; + if(glogp) + glogp->next=this; + prev=glogp; + glogp=this; + } + else{ + if(glogp!=this) // was not set as the global + return -1; + glogp=prev; + if(glogp) + glogp->next=NULL; + prev=NULL; + } + return 1; + } + + virtual void set_throw(int throw_exceptions){ + throw_ex=throw_exceptions; + } + + virtual void set_levels(int out_level=vblALLBAD|vblMESS1, int stop_level=vblFATAL){ + outlevel=out_level; + stoplevel=stop_level; + } + + /// nonzero extra levels are applied instead of set ones + virtual void extra_levels(int out_level=vblALLBAD|vblMESS1, int stop_level=vblFATAL){ + eoutlevel=out_level; + estoplevel=stop_level; + } + + template + int message(const exc_t &signal, int errcode, const char *what, ...){ + int level=log_exception_traits::level(signal); + char buff[1024]; + if(level&(eoutlevel ? eoutlevel : outlevel)){ //needs to print a message + va_list args; + va_start(args,what); + vsnprintf(buff,1024,what,args); + log_text(level,log_exception_traits::name(signal).c_str(),buff); + } + if(level&vblERR){ + if(throw_ex==1) // throws exc_t exception + throw log_exception_traits::add_words(signal,buff); + else if(throw_ex==2) // throws pair<>(int,const char*) exception + throw make_pair(errcode,what); + else if(throw_ex==3) // throws int exception + throw errcode; + } + if(level&(estoplevel ? estoplevel: stoplevel) ){ // needs to stop + exit(errcode); + } + return errcode; + } + + virtual void log_text(int level, const char *messtype, const char *messtext){ + if(descriptor!="") // descriptor is used as header + printf("%s:\n",descriptor.c_str()); + if(messtype!="") + printf("%s: ",messtype); + printf("%s\n",messtext); + } + + /// checks that the deleted one is not in global logger chain + ~message_logger(){ + if(prev){ + prev->next=next; + if(next) + next->prev=prev; + } + set_global(false); + } +}; + +/// global message function +template +int message(const exc_t &signal, int errcode, const char *what, ...){ + if(message_logger::glogp){ + va_list args; + va_start(args,what); + char buff[1024]; + vsnprintf(buff,1024,what,args); + return message_logger::glogp->message(signal,errcode,buff); + } + else + return -1; +} + +/// message logger for which std and error streams may be specified +class stdfile_logger: public message_logger { +protected: + FILE *fout, *ferr; +public: + stdfile_logger(const string &descriptor_="", int throw_exceptions=0, + FILE *out=stdout, FILE *err=stderr, + int out_level=vblALLBAD|vblMESS1,int stop_level=vblFATAL, + int use_globally=0) + : message_logger(descriptor_,out_level,stop_level,throw_exceptions,use_globally),fout(NULL), ferr(NULL){ + set_out(out); + set_err(err); + } + virtual void set_out(FILE *out, int close_prev=0){ + if(close_prev && fout && fout!=stderr && fout !=stdout) + fclose(fout); + fout=out; + } + + virtual void set_err(FILE *err, int close_prev=0){ + if(close_prev && ferr && ferr!=stderr && ferr !=stdout) + fclose(ferr); + ferr=err; + } + + virtual void log_text(int level, const char *messtype, const char *messtext){ + FILE *f= (level&vblALLBAD ? ferr : fout); + if(!f) + return; + if(descriptor!="") // descriptor is used as header + fprintf(f,"%s:\n",descriptor.c_str()); + if(string(messtype)!="") + fprintf(f,"%s: ",messtype); + fprintf(f,"%s\n",messtext); + } +}; + +/// format a string +const char *fmt(const char *format,...); + +/// macros with common usage +#define LOGFATAL(code,text,lineinfo) ((lineinfo) ? ::message(vblFATAL,(code)," %s at %s:%d",(text),__FILE__,__LINE__) : \ + (::message(vblFATAL,(code)," %s",(text))) ) + + +#define LOGERR(code,text,lineinfo) ((lineinfo) ? ::message(vblOERR,(code)," %s at %s:%d",(text),__FILE__,__LINE__) : \ + (::message(vblOERR,(code)," %s",(text))) ) + + +#define LOGMSG(cat,text,lineinfo) ((lineinfo) ? ::message((cat),0," %s at %s:%d",(text),__FILE__,__LINE__) : \ + (::message((cat),0," %s",(text))) ) + + + +# if 0 /// did not work + +/// this may be used to inherit exceptions +/// where level and name are defined whithin a class +struct log_exception { + /// exeption level according to the vbLEVELS + static int level(const log_exception &signal){ return vblFATAL; } + /// the string name of exception category + static string name(const log_exception &signal){ return "undefined exception";} +}; + +/// log_exceptions act as themselves +template<> +struct log_exception_traits{ + static int level(const log_exception &signal){ return log_exception::level(signal); } + static string name(const log_exception &signal){ return log_exception::name(signal); } +}; + +# endif + + +# endif diff --git a/lib/awpmd/ivutils/include/pairhash.h b/lib/awpmd/ivutils/include/pairhash.h new file mode 100644 index 0000000000..f72c544af8 --- /dev/null +++ b/lib/awpmd/ivutils/include/pairhash.h @@ -0,0 +1,683 @@ +/*e*************************************************************************** + * + * Copyright (c), Ilya Valuev 2005 All Rights Reserved. + * + * Author : Ilya Valuev, MIPT, Moscow, Russia + * + * Project : ivutils + * + * + *****************************************************************************/ +/*e**************************************************************************** + * $Log: pairhash.h,v $ + * Revision 1.3 2011/06/11 18:18:50 morozov + * USER-AWPMD compiles on Linux now! + * + * Revision 1.2 2011/06/11 16:53:55 valuev + * sync with LAMMPS + * + * Revision 1.1 2011/06/10 17:15:07 morozov + * First Windows project with the correct directory structure + * + * Revision 1.27 2011/06/09 22:55:08 valuev + * norm matrices + * + * Revision 1.26 2011/06/07 17:43:00 valuev + * added Y derivatives + * + * Revision 1.25 2011/05/24 19:54:32 valuev + * fixed sqmatrix::iterator + * + * Revision 1.24 2011/05/21 23:06:49 valuev + * Norm matrix transform to pysical variables + * + * Revision 1.23 2009/09/24 10:06:38 valuev + * moved matrix printing to template function, reproducing old TB calculations + * + * Revision 1.22 2009/02/10 14:20:45 valuev + * sync with FDTD project + * + * Revision 1.4 2009/01/30 13:54:05 valuev + * restructured as a library + * + * Revision 1.21 2008/08/27 13:34:32 valuev + * made icc-compilable + * + * Revision 1.20 2008/08/25 21:06:11 valuev + * moved using delaration to public + * + * Revision 1.19 2008/07/23 16:55:05 morozov + * *** empty log message *** + * + * Revision 1.18 2008/07/23 16:21:52 morozov + * Corrected Makefile for unix compilation of tcpengine + * + * Revision 1.17 2008/07/18 18:15:31 morozov + * *** empty log message *** + * + * Revision 1.16 2008/07/02 13:11:32 valuev + * new C60+O2 experiments + * + * Revision 1.15 2008/06/24 08:50:00 valuev + * made icc-compilable + * + * Revision 1.14 2008/06/24 08:39:57 valuev + * added ESSL support to TB + * + * Revision 1.13 2008/05/29 14:47:33 valuev + * made icc-compilable + * + * Revision 1.12 2008/05/14 17:17:22 morozov + * Passed 2- and 3-electron test. Added Norm matrix. + * + * Revision 1.11 2008/05/05 17:27:43 morozov + * cvector_3.h is the new header for class cVector_3. Old one is moved to cvector_3old.h + * class hmatrix is added to pairhash.h + * wavepackets.h contains class WavePacket + * + * Revision 1.10 2008/04/21 22:42:30 valuev + * *** empty log message *** + * + * Revision 1.9 2008/04/15 13:11:41 valuev + * Added antisymmetrized wave packets + * + * Revision 1.8 2008/02/28 13:26:04 valuev + * vasp scanner + * + * Revision 1.7 2007/12/13 19:48:59 valuev + * added newlines + * + * Revision 1.6 2006/12/20 14:29:33 valuev + * Updated workflow, sync with FDTD + * + * Revision 1.3 2006/10/27 20:41:01 valuev + * Added detectors sceleton. Updated some of ivutils from MD project. + * + * Revision 1.5 2006/09/26 10:59:42 valuev + * Added nonorthogonal TB (Menon-Subbaswamy) + * + * Revision 1.4 2006/07/21 16:22:03 valuev + * Added Tight Binding for graphite+O + * + * Revision 1.3 2006/04/26 12:12:01 valuev + * Fixed Neighbour Lists (double-single counting), added twostep NL scheme, added Step2 to mdtutorial (use of template potentials), added DelAtom to mdStructure + * + * Revision 1.2 2005/12/09 21:06:38 valuev + * Added neighbour list to mdPotential interface. + * Added missing files to ivutils directory. + * Added mdtutorial and step1 project + * + * Revision 1.1 2005/12/02 18:51:06 valuev + * added HEAD project tree + * + * Revision 1.1 2005/11/30 23:36:11 valuev + * put ivutils to cvs on biolab1.mipt.ru + * + * Revision 1.1 2005/11/30 23:15:43 valuev + * put ivutils on cvs biolab1.mipt.ru + * + * +*******************************************************************************/ +# ifndef PAIRHASH_H +# define PAIRHASH_H + +/*e @file pairhash.h @brief pair hash table +*/ + + +/*r @file pairhash.h @brief работа с хеш-таблицами парных величин +*/ + +# include "refobj.h" + + +///\en Rectangular matrix +template +class recmatrix{ +protected: + mngptr parr; +public: + class iterator{ + friend class recmatrix; + T *ptr; + size_t incr; + iterator(const recmatrix *parent,size_t first_,size_t second_, bool inc_first=false){ + ptr=parent->arr+parent->index(first_,second_); + incr=inc_first ? parent->sizex : 1 ; + }; + iterator(T *ptr_,size_t incr_):ptr(ptr_),incr(incr_){} + public: + iterator(const iterator &other):ptr(other.ptr),incr(other.incr){ + } + iterator():ptr(NULL),incr(0){} + iterator &operator++(){ // prefix + ptr+=incr; + return *this; + } + iterator operator++(int){ // postfix + iterator tmp=*this; + ++*this; + return tmp; + } + iterator operator+(int delta) const { + return iterator(ptr+delta*incr,incr); + } + bool operator!=(const iterator &other) const { + if(ptr!=other.ptr)return true; + else return false; + } + T &operator*() const { + return *ptr; + } + }; + + + T *arr; + size_t sizex, sizey; + + //e default constructor + recmatrix(): parr(NULL,1) { + sizey=sizex=0; + arr=NULL; + } + + //e copy constructor: makes a managed copy + recmatrix(const recmatrix &other):sizex(0),sizey(0),arr(NULL){ + *this=other; + } + + recmatrix &operator=(const recmatrix &other){ + if(this!=&other){ + if(other.sizex*other.sizey<=sizex*sizey) + init(other.sizex,other.sizey,-1); // keeping old array + else + init(other.sizex,other.sizey,1); + size_t n=get_datasize(sizex,sizey); + for(size_t i=0;i=0){ // for changing the managed flag? + parr.reset(parr.ptr(),smanaged ? smanaged|0x8 : 0 ); + managed=smanaged; + } + if(sizex==nx && sizey==ny) // no need to allocate + return 1; + sizex=nx; + sizey=ny; + if(managed){ + if(sizex>0 && sizey>0) + parr.reset(new T[get_datasize(sizex,sizey)],managed|0x8); + arr=parr.ptr(); + } + return 1; + } + + recmatrix(size_t nx, size_t ny):sizex(0), sizey(0){ + init(nx,ny,1); + } + + //e initializes by unmanaged pointer + recmatrix(size_t nx, size_t ny , T *ptr):parr(ptr,0),sizex(nx), sizey(ny) { + init(nx,ny); + } + + //e attaches to new pointer and sets unmanaged size + void AttachTo(size_t nx,size_t ny, T *ptr){ + init(0,0); + sizex=nx; + sizey=ny; + parr.reset(ptr,0); + } + + void Set(const T &val){ + size_t i, n=get_datasize(sizex,sizey); + for(i=0;isizey? sizey: sizex); + for(i=0;i +class sqmatrix: public recmatrix { + +public: + + size_t size; + + //e default constructor + sqmatrix(){} + + //e copy constructor: makes a managed copy + sqmatrix(const sqmatrix &other):size(0){ + *this=other; + } + + sqmatrix &operator=(const sqmatrix &other){ + if(this!=&other){ + *((recmatrix *)this)=*((recmatrix *)&other); + size=other.size; + } + return *this; + } + + virtual size_t get_datasize(size_t n) const{ + return n*n; + } + + + virtual int init(size_t n, int smanaged=-1){ + size=n; + return recmatrix::init(n,n,smanaged); + int managed=recmatrix::parr.managed(); + } + + sqmatrix(size_t n):size(0){ + init(n,1); + } + + //e initializes by unmanaged pointer + sqmatrix(size_t n, T *ptr):size(n){ + init(n); + } + + //e attaches to new pointer and sets unmanaged size + void AttachTo(size_t n, T *ptr){ + init(0); + size=n; + recmatrix::parr.reset(ptr,0); + } + + +}; + + + +# if 0 +//e square matrix +template +class sqmatrix{ + mngptr parr; +public: + class iterator{ + friend class sqmatrix; + T *ptr; + size_t incr; + iterator(const sqmatrix *parent,size_t first_,size_t second_, bool inc_first=false){ + ptr=parent->arr+parent->index(first_,second_); + incr=inc_first ? parent->size : 1 ; + }; + iterator(T *ptr_,size_t incr_):ptr(ptr_),incr(incr_){} + public: + iterator(const iterator &other):ptr(other.ptr),incr(other.incr){ + } + iterator():ptr(NULL),incr(0){} + iterator &operator++(){ // prefix + ptr+=incr; + return *this; + } + iterator operator++(int){ // postfix + iterator tmp=*this; + ++*this; + return tmp; + } + iterator operator+(int delta) const { + return iterator(ptr+delta*incr,incr); + } + bool operator!=(const iterator &other) const { + if(ptr!=other.ptr)return true; + else return false; + } + T &operator*() const { + return *ptr; + } + }; + + + T *arr; + size_t size; + + //e default constructor + sqmatrix(): parr(NULL,1) { + size=0; + arr=NULL; + } + + //e copy constructor: makes a managed copy + sqmatrix(const sqmatrix &other):size(0),arr(NULL){ + *this=other; + } + + sqmatrix &operator=(const sqmatrix &other){ + if(this!=&other){ + if(other.size<=size) + init(other.size,-1); // keeping old array + else + init(other.size,1); + size_t n=get_datasize(size); + for(size_t i=0;i=0){ // for changing the managed flag? + parr.reset(parr.ptr(),smanaged ? smanaged|0x8 : 0 ); + managed=smanaged; + } + if(size==n) // no need to allocate + return 1; + size=n; + if(managed){ + if(size>0) + parr.reset(new T[get_datasize(size)],managed|0x8); + arr=parr.ptr(); + } + return 1; + } + + sqmatrix(size_t n):size(0){ + init(n,1); + } + + //e initializes by unmanaged pointer + sqmatrix(size_t n, T *ptr):parr(ptr,0),size(n){ + init(n); + } + + //e attaches to new pointer and sets unmanaged size + void AttachTo(size_t n, T *ptr){ + init(0); + size=n; + parr.reset(ptr,0); + } + + void Set(const T &val){ + size_t i, n=get_datasize(size); + for(i=0;i +int fileout(FILE *f, const matrix_t &matr, const char *elm_fmt, const char *elm_sep=" ", const char *line_sep="\n"){ + size_t i, j; + int res=0; + for(i=0;i +class smatrix: public sqmatrix{ + typedef sqmatrix base_t; +public: + + virtual size_t get_datasize(size_t n) const{ + return n*(n+1)/2; + } + + size_t index(size_t i, size_t j) const { + if(i>=j) + return (2*base_t::size-j-1)*j/2+i; + else + return (2*base_t::size-i-1)*i/2+j; + } + + T &operator()(size_t i,size_t j){ + return base_t::arr[index(i,j)]; + } + + T operator()(size_t i,size_t j) const { + return base_t::arr[index(i,j)]; + } + + void set(long i,long j,const T& val){ + (*this)(i,j)=val; + } + + void SetDiag(const T &val){ + size_t i; + for(i=0;i(other){} + + smatrix &operator=(const smatrix &other){ + return (smatrix&)( *(sqmatrix *)this = (sqmatrix &)other ); + } + + smatrix(size_t n):sqmatrix(n){} + + //e initializes by unmanaged pointer + smatrix(size_t n, T *ptr):sqmatrix(n,ptr){} + +}; + + + +//e Hermitian matrix +template +class hmatrix : public smatrix +{ +public: + using smatrix::arr; + using smatrix::size; + hmatrix() : smatrix() {} + + hmatrix(const smatrix &other) : smatrix(other) {} + + //e copy constructor: makes a managed copy + hmatrix(const hmatrix &other): smatrix(other){} + + hmatrix &operator=(const hmatrix &other) { + return (hmatrix&)( *(smatrix*)this = (smatrix&)other ); + } + + hmatrix(size_t n) : smatrix(n) {} + + hmatrix(size_t n, T *ptr) : smatrix(n, ptr) {} + + T operator()(size_t i, size_t j) const { + if(i<=j) return arr[(2*size-i-1)*i/2+j]; + else return conj( arr[(2*size-j-1)*j/2+i] ); + } + + void set(long i,long j,const T& val){ + if(i<=j) arr[(2*size-i-1)*i/2+j] = val; + else arr[(2*size-j-1)*j/2+i] = conj(val); + } + +}; + + +//e Basic pair hash class +template +class PairHash{ +public: + //e find the value with indexes i, j + //e @return 0 if not found, 1 otherwise + //e if retval is not NULL, puts the found value there + virtual int Find(long i, long j, T *retval=NULL)=0; + virtual int Find(long i, long j, T **retval=NULL)=0; + virtual int Del(long i, long j)=0; + virtual int Put(long i, long j, const T *value)=0; + virtual int Put(long i, long j, const T& value)=0; + virtual int Clear()=0; + virtual ~PairHash(){} +}; + +//e Hash with symmetric matrix +template +class PairHashM: public PairHash{ + smatrix indm; + T *arr; + int as; +public: + PairHashM(long n, int antisymmetric =0):indm(n),as(antisymmetric){ + indm.Set(-1); + arr= new T[n*(n+1)/2]; + } + int Find(long i, long j, T *retval=NULL){ + long ind=indm(i,j); + if(ind>=0){ + if(retval){ + if(as && i=0){ + *retval=&arr[ind]; + if(as && i &arg) is added + * + * Revision 1.31 2009/01/30 13:54:05 valuev + * restructured as a library + * + * Revision 1.30 2009/01/21 09:28:15 lesha + * refvector::clear is added + * + * Revision 1.29 2009/01/15 07:31:07 lesha + * *** empty log message *** + * + * Revision 1.28 2009/01/14 10:02:36 lesha + * operator [] is added to mngptr + * + * Revision 1.27 2008/04/29 01:23:59 lesha + * nothing important + * + * Revision 1.26 2008/02/28 08:57:27 lesha + * shptr::free() is made public + * + * Revision 1.25 2008/02/27 13:37:23 lesha + * shptr is added + * + * Revision 1.24 2008/01/22 10:14:05 lesha + * mngarg is added + * + * Revision 1.23 2007/08/08 10:55:37 lesha + * constructor in refvector is correted + * + * Revision 1.22 2007/07/10 19:52:44 lesha + * make gcc compilable + * + * Revision 1.21 2007/07/06 12:23:37 valuev + * made compilable with icc 9 + * + * Revision 1.20 2007/06/22 09:42:25 valuev + * *** empty log message *** + * + * Revision 1.19 2007/06/17 00:51:44 lesha + * refobj, refcounter :: reset is modified for ptr==this->ptr case + * + * Revision 1.18 2007/06/05 16:30:53 lesha + * make gcc compilable + * + * Revision 1.17 2007/06/05 11:37:53 lesha + * make gcc compilable + * + * Revision 1.16 2007/06/05 11:07:04 lesha + * make gcc compilable + * + * Revision 1.15 2007/06/04 14:03:55 lesha + * *** empty log message *** + * + * Revision 1.14 2007/05/31 18:00:42 lesha + * *** empty log message *** + * + * Revision 1.13 2007/05/31 16:57:30 lesha + * ref_sequence is added + * + * Revision 1.12 2007/05/31 01:25:01 lesha + * new version of mng_ptr, pencil etc + * + * Revision 1.11 2007/02/20 10:26:11 valuev + * added newlines at end of file + * + * Revision 1.10 2007/02/16 10:16:51 valuev + * allowed array mngptr + * + * Revision 1.4 2007/02/16 09:40:32 valuev + * Added Nudged Elastic Band saddle point search + * + * Revision 1.3 2006/12/20 14:29:33 valuev + * Updated workflow, sync with FDTD + * + * Revision 1.9 2006/12/14 08:42:36 valuev + * reformulated detector + * projectors, corrected open file limit control, tested Fourier sceleton + * + * Revision 1.8 2006/11/29 18:05:05 valuev + * made the code compilable with g++ + * + * Revision 1.7 2006/11/29 17:17:01 valuev + * added using base_t::member for ANSI-compatibility + * + * Revision 1.6 2006/11/29 17:11:59 valuev + * added includes + * + * Revision 1.5 2006/11/28 09:16:59 valuev + * Fixed vectors storing managed pointers + * + * Revision 1.4 2006/11/24 20:17:31 valuev + * Added CVS headers + * +*******************************************************************************/ +#ifndef _REFOBJ_H +#define _REFOBJ_H + +# include +# include +# include + +using namespace std; + +template +class mngarg: public pair{ +public: + typedef pair base_t; + using base_t::second; + using base_t::first; + + mngarg(T *ptr, int managed=0): pair(ptr,managed){} + template + mngarg(const mngarg &arg): pair(arg.first,arg.second){} +}; + +template +mngarg make_mngarg(T *ptr, int managed=1){ + return mngarg(ptr,managed); +} + +/// managed pointer +/// managed==0 do not delete +/// managed==1 delete +/// (NOT IMPLEMENTED) managed==2 copy and delete, requires copy constructor +/// flag 0x8 -- delete as array +template +class mngptr: public pair{ +public: + typedef pair base_t; + typedef T *pointer; + using base_t::second; + using base_t::first; + + mngptr(T* ptr=NULL, int managed=0): pair(ptr,managed){ + //if(managed==2)ptr= new T(*ptr); + } + mngptr(const mngarg &arg): pair(arg.first,arg.second){} + const mngptr &operator=(const mngarg &arg){ + reset(arg.first,arg.second); + return *this; + } + void reset(T* ptr=NULL, int managed=0){ + if(second && first && first!=ptr){ + if(second&0x8)delete [] first; + else delete first; + } + first=ptr; + second=managed; + } + void reset(const mngarg &arg){ + reset(arg.first,arg.second); + } + T* ptr() const { + return first; + } + T* operator->() const { + return first; + } + T& operator*() const{ + return *first; + } + T& operator[] (int i) const{ + return *(first+i); + } + int managed() const { + return second; + } + ~mngptr(){ + reset(); + } +}; + + +# if 0 +template class cont_tt, class T> +class refcontainer: public cont_tt< T * >{ +protected: + int man; +public: + typedef cont_tt< T * > base_t; + typedef typename base_t::iterator iterator; + typedef typename base_t::const_iterator const_iterator; + + refcontainer(int smanaged=0):man(smanaged){} + refcontainer(size_t n, int smanaged=0):base_t(n),man(smanaged){} + + void set_managed(int sman){ + man=sman; + } + + ~refcontainer(){ + if(man){ + size_t i, n=base_t::size(); + for(i=0;i +class refvector: public std::vector< T * >{ +protected: + int man; +public: + typedef vector base_t; + typedef typename base_t::iterator iterator; + typedef typename base_t::const_iterator const_iterator; + + refvector(int smanaged=0):man(smanaged){} +// refvector(size_t n, int smanaged=0):base_t(n),man(smanaged){} // ambigious constructors + refvector(size_t n, int smanaged):base_t(n),man(smanaged){} + + void set_managed(int sman){ + man=sman; + } + + ~refvector(){ + clear(); + } + + void clear() { + if(man){ + iterator it=base_t::begin(); + for(;it!=base_t::end();++it) + if(*it) + delete (*it); + } + base_t::clear(); + } + + iterator erase(iterator it){ + if(man && *it) + delete (*it); + return base_t::erase(it); + } + +}; + + +template +class refmap: public std::map{ +protected: + int man; +public: + typedef std::map base_t; + typedef typename base_t::iterator iterator; + typedef typename base_t::const_iterator const_iterator; + + refmap(int smanaged=0):man(smanaged){} + refmap(size_t n, int smanaged=0):base_t(n),man(smanaged){} + + void set_managed(int sman){ + man=sman; + } + + ~refmap(){ + clear(); + } + + void clear() { + if(man){ + for(typename base_t::iterator i=base_t::begin();i!=base_t::end();++i) + if(i->second) delete i->second; + } + base_t::clear(); + } + + iterator erase(iterator it){ + if(man && it->second) + delete it->second; + return base_t::erase(it); + } +}; + + +template +class delete_ptr{ +public: + void operator()(T *ptr){ + delete ptr; + } +}; + +template > +class shptr{ + template friend class shptr; + T *p; + int *num; //if num==NULL than p is not managed (as in mngptr) + + void set(T *p_, int managed){ + p=p_; + if(p&&managed){ + num=new int; + *num=1; + } + else num=NULL; + } + template + void set(const Y &other){ + p=other.p; + if(p){ + num=other.num; + if(num)(*num)++; + } + else num=NULL; + } + void set(const shptr &other){ + p=other.p; + if(p){ + num=other.num; + if(num)(*num)++; + } + else num=NULL; + } + +public: + shptr(T* p=NULL, int managed=1){ + set(p,managed); + } + shptr(const mngarg &arg){ + set(arg.first,arg.second); + } + template + shptr(const Y &other){ + set(other); + } + shptr(const shptr &other){ + set(other); + } + + void reset(T *p_, int managed=1) { + if(p!=p_){ + free(); + set(p_,managed); + } + } + void reset(const shptr &other) { + if(this!=&other){ + free(); + set(other); + } + } + + const shptr &operator=(T *p) { + reset(p,0); + return *this; + } + const shptr &operator=(const mngarg &arg) { + reset(arg.first,arg.second); + return *this; + } + template + const shptr &operator=(const Y &other){ + reset(other); + return *this; + } + const shptr &operator=(const shptr &other) { + reset(other); + return *this; + } + + virtual ~shptr(){ + free(); + } + void free(){ + if(p){ + if(num){ + (*num)--; + if((*num)==0){ + delete_t()(p); + delete num; + } + num=NULL; + } + p=NULL; + } + } + + bool valid() const { + return p!=NULL; + } + + T* ptr() const { + return p; + } + T* operator->() const { + return p; + } + T& operator*() const{ + return *p; + } +}; + + + + +/* +class RefObject{ + void *ref_data; + int ref_count; +public: + +protected: + virtual void delete_data(void *data); + virtual void *new_data(); + virtual void *copy_data(void *data); +} + + + +class RefA: public RefObject{ + +public: + refA(){ + ref_data = new A; + } + refA(const refA &other){ + Ref(other); + } + refA& operator=(const refA &other){ + if(ref_data != other.ref_data){ + Ref(other); + } + return *this; + } +private: + void delete_data(void *data){ + delete (A *)data; + } + void *new_data(){ + return (void *)new A; + } + void *copy_data(void *data){ + return (void *)(new A(*((A*)data))); + } + + + +}*/ + + +#endif diff --git a/lib/awpmd/ivutils/include/vector_3.h b/lib/awpmd/ivutils/include/vector_3.h new file mode 100644 index 0000000000..c269628c32 --- /dev/null +++ b/lib/awpmd/ivutils/include/vector_3.h @@ -0,0 +1,707 @@ +/*s*************************************************************************** + * + * Copyright (c), Ilya Valuev 2005 All Rights Reserved. + * + * Author : Ilya Valuev, MIPT, Moscow, Russia + * + * Project : GridMD, ivutils + * + *****************************************************************************/ + +/*s**************************************************************************** + * $Log: vector_3.h,v $ + * Revision 1.2 2011/06/11 16:53:55 valuev + * sync with LAMMPS + * + * Revision 1.1 2011/06/10 17:15:07 morozov + * First Windows project with the correct directory structure + * + * Revision 1.22 2010/11/17 02:13:32 valuev + * Added analysis phase, fixed some memory leaks + * + * Revision 1.21 2009/09/09 14:43:35 valuev + * fixed trajreader + * + * Revision 1.20 2009/07/24 05:08:46 valuev + * Sync with FDTD, added molecule setup + * + * Revision 1.33 2009/06/01 13:01:42 valuev + * Added ShearBox + * + * Revision 1.32 2009/05/28 07:49:00 valuev + * updated chemosensor + * + * Revision 1.31 2009/05/22 08:53:52 valuev + * new uiExperiment interface + * + * Revision 1.30 2009/02/22 10:49:56 lesha + * Vector_Nt constructors are modified + * + * Revision 1.29 2009/02/04 14:23:31 valuev + * fixed bug in maxabscoord and minabscoord functions + * + * Revision 1.28 2009/02/04 10:56:30 lesha + * SINGLE_PRECISION is recovered + * + * Revision 1.27 2009/02/03 00:47:35 lesha + * *** empty log message *** + * + * Revision 1.26 2009/01/30 13:54:05 valuev + * restructured as a library + * + * Revision 1.16 2008/08/18 21:40:09 valuev + * added Gurski-Krasko potential + * + * Revision 1.15 2008/05/05 17:27:43 morozov + * cvector_3.h is the new header for class cVector_3. Old one is moved to cvector_3old.h + * class hmatrix is added to pairhash.h + * wavepackets.h contains class WavePacket + * + * Revision 1.14 2008/04/22 12:44:17 valuev + * made gcc 4.12 compilable + * + * Revision 1.13 2008/04/21 23:13:44 valuev + * made gcc 4.12 compilable + * + * Revision 1.12 2008/04/21 22:42:30 valuev + * *** empty log message *** + * + * Revision 1.11 2008/04/15 13:11:41 valuev + * Added antisymmetrized wave packets + * + + * +*******************************************************************************/ +/*r @file vector_3.h @brief работа с трехмерными векторами +*/ + +# ifndef VECTOR_3_H + +# define VECTOR_3_H + +# define _USE_MATH_DEFINES + +# include +# include +# include +# include + +// some compilers don't define PI! +# ifndef M_PI +# define M_PI 3.1415926535897932385 +# endif + +# ifndef fmod +//r деление по модулю чисел с плавающей точкой +# define fmod(a,b) ((a)-((long)((a)/(b))*(b))) +# endif + +using namespace std; + +#ifndef SINGLE_PRECISION +typedef double vec_type; +#else +typedef float vec_type; +#endif + +//e "infinitely" large number in Vector_3 sense +//r "бесконечно большое" число в смысле Vector_3 +//# ifndef SINGLE_PRECISION +//# define VEC_INFTY 1.e20 +//# else +//# define VEC_INFTY 1.e20f +//# endif +#define VEC_INFTY numeric_limits::max() + +//e "infinitely" small number in Vector_3 sense (used for vector comparisons) +//r "бесконечно малое" число в смысле Vector_3 (используется для сравнений) +//# ifndef SINGLE_PRECISION +//# define VEC_ZERO 1.e-20 +//# else +//# define VEC_ZERO 1.e-20f +//# endif +#define VEC_ZERO 512*numeric_limits::epsilon() + +//r N-мерный вектор типа T, с некоторыми полезными операциями +template +struct Vector_Nt { + typedef T value_type; + + T v[N]; + + Vector_Nt(const T &a=0){ + for (int i=0; i0)v[0]=a1; + if(N>1)v[1]=a2; + for(int i=2;i0)v[0]=s1; + if(N>1)v[1]=s2; + if(N>2)v[2]=s3; + for(int i=3;i + Vector_Nt(const A *beg) { + for (int i=0; i + Vector_Nt(const Vector_Nt& v1) { + for (int i=0; i + void copy_to(A *beg) const { + for (int i=0; iVEC_ZERO)return false; + return true; + }; + + inline bool operator!=(const Vector_Nt &vect) const{ + return (!(this->operator==(vect))); + }; + + inline Vector_Nt operator+(const Vector_Nt& vect) const{ + Vector_Nt result; + for (int i=0; inorm(); + if(norm>=VEC_ZERO){ + T c=newnorm/norm; + for (int i=0; i=VEC_INFTY){ + if(result>=0) + result=0.; + result-=1; + v[i]=v[i]>0 ? 1 : -1; + } + else if(result>=0) //still summing the normal components + result+=v[i]*v[i]; + else + v[i]=0.; + } + if(fabs(result)cell[i]/2){ + ret[i]-=cell[i]; + } + else if(v[i]<-cell[i]/2){ + ret[i]+=cell[i]; + } + } + } + return ret; + } + + //e @return a vector projection on a given axis + Vector_Nt prj(int i) const { + Vector_Nt res; + res[i]=v[i]; + return res; + } + + //e @return a vector projection on a given vector (k*v)*k + Vector_Nt prj(const Vector_Nt &k) const { + return k*(k*(*this)); + } + + //e @return a vector of length l parallel to this + Vector_Nt unit(T l=1) const { + Vector_Nt res(*this); + res.normalize(l); + return res; + } + + //e reduction to elementary cell [0, cell[i]) (FOR REDUCTION TO ELEMENTARY CELL) + //e flags indicate the periodicity in specific directions: 0x1 for X, 0x2 for Y, 0x4 for Z + //r Почти то же, что и rcell1, но без ограничения на положение *this и с другой системой ячеек. + /*r В начале координат находится не центр ячейки, а ее угол. Может работать медленнее из-за наличия + операции деления по модулю с плавающей точкой */ + Vector_Nt rcell(const Vector_Nt &cell, int flags=0xffff) const { + Vector_Nt ret(*this); + for (int i=0, flag=1; icell[i]) + // printf("!"); + } + } + return ret; + } + + Vector_Nt rpcell(const Vector_Nt &p1, const Vector_Nt &cell, int flags=0xfff) const { + Vector_Nt ret(*this); + for (int i=0, flag=1; ip1[i]+cell[i]) { + ret.v[i]=fmod(v[i]-p1[i],cell[i])+p1[i]; + if (ret.v[i]vv){ + im=i; + vv=v[i]; + } + } + if(ind)*ind=im; + return vv; + } + + + //e returns the corrd having maximal absolute value + T maxabscoord(int *ind=NULL) const { + int im=0; + T vv=fabs(v[0]); + for (int i=1; ivv){ + im=i; + vv=fabs(v[i]); + } + } + if(ind)*ind=im; + return v[im]; + } + + //e returns the corrd having minimal absolute value + T minabscoord(int *ind=NULL) const { + int im=0; + T vv=fabs(v[0]); + for (int i=1; i=VEC_INFTY) + return true; + } + return false; + } +}; + +template +Vector_Nt operator*(const T &coeff,const Vector_Nt &vec){ + return vec*coeff; +} + +template +Vector_Nt operator*(int coeff,const Vector_Nt &vec){ + return vec*coeff; +} + +//template <> Vector_Nt operator*(const double &coeff,const Vector_Nt &vec); + +// old Vector_3 compatibility typedefs and functions +typedef Vector_Nt Vector_3; +typedef Vector_3 *Vector_3P; +typedef Vector_Nt Vector_2; +typedef Vector_2 *Vector_2P; + +template +class Vector_N: public Vector_Nt{ +}; + +//Vector_3 operator*(int coeff,const Vector_3 &vec){ +// return vec*((vec_type)coeff); +//} + + +//e finds the maximum distance between vector pairs +//r Находит максимальное расстояние между векторами va1[i], va2[i], i=1..n +/*r @param va1 - массив Vector_3[n] + @param n - длина массивов va1 и va2 +*/ +vec_type dist_max(Vector_3 *va1,Vector_3 *va2,int n); + +//e finds average distance between vector pairs +//r Находит среднее расстояние между векторами va1[i], va2[i], i=1..n +vec_type dist_av(Vector_3 *va1,Vector_3 *va2,int n); + +//e finds the average difference norm between two vector sets of the same length +/*e optionally gives the indexes for maximal and minimal difference + va2 can be NULL, then the norm of va1 is used */ + +//r Находит среднее расстояние между va1[i] и va2[i], а также, по желанию, индексы, на которых достигается min и max расстояние +vec_type diff_av(Vector_3 *va1,Vector_3 *va2,int n, int *minind=0, int *maxind=0); + +//e finds suitable perpendicular to a vector +//r Находит перпендикуляр к вектору vAB +Vector_3 FindPerp(const Vector_3 &vAB); + + + +//e Returns the average (center) vector of the vector array +//e and cooordinates of a minimal box in which it is contained +Vector_3 GetScope(const Vector_3 *varr,long n,Vector_3* box_min,Vector_3* box_max); + +//e the same with long index array +Vector_3 GetIScope(const Vector_3 *varr,long *indarr,long n,Vector_3* box_min,Vector_3* box_max); + +//e the same with int index array +Vector_3 GetIScopei(const Vector_3 *varr,int *indarr,int n,Vector_3* box_min,Vector_3* box_max); + +// neue Funktionen + +//e clears vector array with optional integer index +//r Очистка массива векторов, с поддержкой индексирования +/*r +В данном Vector_3 vec[] обнуляет n координат. Если ind==NULL, то +очищает первые n элементов. Если ind!=NULL, то для i=0..n-1 +очищает vec[ind[i]] +См. @ref indexed_calculations. +*/ +void clear_vecarri(int n,Vector_3 *vec, int *ind=0); + +//e reflects the vector ini+dir*t+0.5*force*t^2 to be inside a box limited by 0 and box sizes +//e changes dir according to the final state +//e fills crossed dir with bit flags corresponding directions along which the walls were crossed +Vector_3 Reflect(Vector_3& ini, double t,Vector_3 &dir, double *box, int flag=0x7, const Vector_3 &force=Vector_3()); + + +inline vec_type vec_area(const Vector_2 &vect1, const Vector_2 &vect2) { + return fabs(vect1[0]*vect2[1]-vect1[1]*vect2[0])/2; +}; + +inline vec_type vec_area(const Vector_3 &vect1, const Vector_3 &vect2) { + return (vect1%vect2).norm()/2; +}; + +// remake for vec_type! + +template +denum_t acccomp(const num_t x, const denum_t y, const val_t val) { + num_t eps_num=numeric_limits::epsilon(); + denum_t eps_denum=numeric_limits::epsilon(); + return (eps_num>=eps_denum) ? acccomp(x, (num_t)y, (num_t)val) : acccomp((denum_t)x, y, (denum_t)val); +} + +template +denum_t acccomp(const num_t x, const denum_t y) { + return acccomp(x, y, num_t(1)); +} + +template +bool acccomp(const T x, const T y, const T val) { + T eps=numeric_limits::epsilon(); + int iexp; + frexp(val,&iexp); + T mult=ldexp(.5, iexp+1); + T err=T(256)*mult*eps; + return fabs(x-y)<=err; +} + +template +bool acccomp(const T x, const T y) { + return acccomp(x, y, T(1)); +} + +template +denum_t accdiv(const num_t x, const denum_t y) { + num_t eps_num=numeric_limits::epsilon(); + denum_t eps_denum=numeric_limits::epsilon(); + return (eps_num>=eps_denum) ? accdiv1(x, (num_t)y) : accdiv1((denum_t)x, y); +} + + +template +T accdiv1(T x, T y) { + T result; + T eps=numeric_limits::epsilon(); + + T fr=x/y; + T ifr=floor(fr+T(.5)); +// T err=64*eps*(ifr!=0 ? fabs(ifr) : 1); + int mult; + if (ifr<=512) + mult=512; + else { + int iexp; + frexp(ifr,&iexp); + mult=int(ldexp(.5, iexp+1)); + } + T err=mult*eps; + if (fabs(fr-ifr)<=err) + result=ifr; + else + result=fr; + return result; +} + + + + +template +denum_t accdiv_rmd(const num_t x, const denum_t y, P *remainder) { + num_t eps_num=numeric_limits::epsilon(); + denum_t eps_denum=numeric_limits::epsilon(); + return (eps_num>=eps_denum) ? accdiv_rmd(x, (num_t)y, remainder) : accdiv_rmd((denum_t)x, y, remainder); +} + +template +T accdiv_rmd(const T x, const T y, P *remainder) { + T result=accdiv(x, y); + T iresult=floor(result); + if (result-iresult>0) + *remainder=x-iresult*y; + else + *remainder=0; + return result; +} + +//e returns random unit vector uniformely distributed in space (?? check this) +inline Vector_3 randdir(){ + vec_type xi1=2*((vec_type)rand())/RAND_MAX-1.; + vec_type xi2=((vec_type)rand())/RAND_MAX; + vec_type r1=sqrt(1.-xi1*xi1); + return Vector_3(r1*cos(2*M_PI*xi2),r1*sin(2*M_PI*xi2),xi1); +} + +///\en Calculates extent of the vector container. +/// \return the center of the vector set, optionally +/// (if arguments are not NULL) fills the bounding box in \a box_min, \a box_max. +template +Vector_3 get_extent(vec_inp_it beg,vec_inp_it end, Vector_3* box_min=NULL,Vector_3* box_max=NULL){ + if(beg==end) + return Vector_3(); + Vector_3 center(*beg++); + Vector_3 cube1(center), cube2(center); + size_t n=1; + for(;beg!=end;++beg){ + Vector_3 vec=*beg; + center+=vec; + for(size_t j=0;j<3;j++){ + if(cube1[j]>vec[j]) + cube1[j]=vec[j]; + if(cube2[j] +bool much_less(const T& a, const T& b, const T &factor){ + if(fabs(a) +# include +# include +# include "cvector_3.h" + +using namespace std; + +/** @file wpmd.h + @brief Classes to handle Gaussian Wave Packets. */ + +// Constants +const double MIN_EXP_ARG = -15.; // Minimum noticeable argument for exp function + +class WavePacket; + + +///\en Template for v=der operation in \ref Wavepacket::int2phys_der() +template +struct eq_second : public binary_function { + Type operator()(const Type& _Left, const Type& _Right) const{ + return _Right; + } +}; + +///\en Template for v=-der operation in \ref Wavepacket::int2phys_der() +template +struct eq_minus_second : public binary_function { + Type operator()(const Type& _Left, const Type& _Right) const{ + return -_Right; + } +}; + +///\en Compares complex numbers on a per component basis. +/// \return \retval 0 if all component differences are 0 within tolerance \a tol (EQUAL), +/// \retval -1 for LESS +/// \retval 2 for GREATER +template< class CT > +int compare_compl(const CT &a, const CT& b, double tol=0.){ + double dd=real(a)-real(b); + if(dd<-tol) + return -1; + if(dd>tol) + return 1; + dd=imag(a)-imag(b); + if(dd<-tol) + return -1; + if(dd>tol) + return 1; + return 0; +} + + +///\en Compares vectors on a per component basis. +/// \return \retval 0 if all component differences are 0 within tolerance \a tol (EQUAL), +/// \retval -1 for LESS +/// \retval 2 for GREATER +inline int compare_vec(const Vector_3 &a, const Vector_3& b, double tol=0.){ + for(int i=0;i<3;i++){ + double dd=a[i]-b[i]; + if(dd<-tol) + return -1; + if(dd>tol) + return 1; + } + return 0; +} + + +/// wavepacket is w(x)=exp(-a*x^2+b*x+lz) +class WavePacket{ + /// constructs a conjugate packet + friend WavePacket conj(const WavePacket &wp); + +public: + cdouble a; + cVector_3 b; + cdouble lz; + + WavePacket(const cdouble &sa=cdouble(1.,0.),const cVector_3 &sb=cVector_3(0.,0.), const cdouble &slz=0.): a(sa), b(sb), lz(slz){ + } + + WavePacket operator*(const WavePacket& other) const { + return WavePacket(a+other.a,b+other.b,lz+other.lz); + } + + /// returns the integral of w(x) over 3D space + cdouble integral() const { + cdouble z = lz + b.norm2()/(4.*a); + if(real(z) < MIN_EXP_ARG) return 0.; + return pow(M_PI/a,3./2.)*exp(z); + } + + /// init normalized packet with physical parameters: r0, p0, width, pw + /// w(x)=(3/2pi width^(3/4)exp[-(3/(4 width^2)-i pw/(2*width) )(x-r)^2+i p (x-r)] + void init(const double width=1., const Vector_3 &r=0., const Vector_3 &p=0., const double pw=0.){ + a = (3./(2.*width) - cdouble(0.,1.)*pw)/(2.*width); + b = (2.*a)*r + cdouble(0.,1.)*p; + lz = log(3./(2.*M_PI*width*width))*(3./4.) + (-a*r.norm2() - cdouble(0.,1.)*(p*r)); + } + + /// init normalized packet with complex parameters a and b + /// w(x)=(3/2pi width^(3/4)exp[-(3/(4 width^2)-i pw/(2*width) )(x-r)^2+i p (x-r)] + void init(const cdouble &a_, const cVector_3 &b_){ + a = a_; + b = b_; + Vector_3 r = get_r(); + double r2 = r.norm2(); + lz = cdouble( log( real(a)*(2./M_PI) )*(3./4.) - r2*real(a), r2*imag(a) - r*imag(b) ); + } + + cdouble operator()(const Vector_3 &x) const { + return exp(lz - a*x.norm2() + b*x); + } + + /// ajusts lz so that Integral[w*(conj(w))] is 1 after this operation + WavePacket& normalize(){ + //lz=0.; + //lz=overlap(conj(*this)); + //lz=(-1./2)*log(lz); + Vector_3 r = get_r(); + double r2 = r.norm2(); + lz = cdouble( log( real(a)*(2./M_PI) )*(3./4.) - r2*real(a), r2*imag(a) - r*imag(b) ); + return *this; + } + + /// computes 3D overlap of wavepackets Inegral[w*other] + cdouble overlap(const WavePacket &other) const{ + WavePacket wp=(*this)*other; + return wp.integral(); + } + + /// returns translated packet to the position of r'=r+dr + WavePacket translate(const Vector_3 &dr) const { + WavePacket wp(a,b,lz); + wp.b+=2*dr*a; + Vector_3 r=get_r(); + double dr2=(r+dr).norm2()-r.norm2(); + wp.lz+=-a*dr2-cdouble(0.,(get_p()*dr)); + return wp; + } + + /// width + double get_width() const{ + return sqrt(3./4./real(a)); + } + /// width momentum + double get_pwidth() const{ + return -imag(a)*sqrt(3./real(a)); + } + /// both width and width momentum + pair get_width_pars() const{ + double c=sqrt(3./2./real(a)); + return make_pair(c/2,-imag(a)*c); + } + /// position + Vector_3 get_r() const { + return real(b)/(2*real(a)); + } + /// momentum + Vector_3 get_p() const { + return imag(b) - real(b)*(imag(a)/real(a)); + } + + ///\en Transforms derivatives of a function whith respect to WP parameters + /// from internal into physical representation, i. e.:\n + /// from df/d{are,aim,b0re,b0im,b1re,b1im,b2re,b2im} (8 values accessed by input iterator d_it in the given order)\n + /// to df/d{x0,x1,x2}, df/d{p0,p1,p2}, df/dw, df/dpw + /// The supplied inputs (val) are modified by op: val=op(val,phys_der). + /// Use operation=eq_second for the supplied inputs to be replaced by new physical derivative values. + /// The inpput and output locations may coinside, an internal buffer is used for transformation. + template class operation, class d_it, class dfdx_it, class dfdp_it, class dfdw_it, class dfdpw_it> + void int2phys_der(d_it dfdi_,dfdx_it dfdx, dfdp_it dfdp, dfdw_it dfdw, dfdpw_it dfdpw, double h_p=h_plank) const { + operation op; + double dfdi[8], dfn[8];// internal buffer + for(int i=0;i<8;i++) + dfdi[i]=*dfdi_++; + double w=get_width(); + Vector_3 r=get_r(); + double t=3/(2*w*w*w); + dfn[6]= -t*dfdi[0]-imag(a)*dfdi[1]/w; //dw + dfn[7]=-dfdi[1]/(2*w*h_p); // dpw + for(int i=0;i<3;i++){ + dfn[i]= 2*real(a)*dfdi[2+2*i]+2*imag(a)*dfdi[2+2*i+1]; + dfn[3+i]= dfdi[2+2*i+1]*(/*m_electron*/1./h_p) ; //*(h_plank/m_electron); + dfn[7]+=-(r[i]*dfdi[2+2*i+1]/w)/h_p; + dfn[6]+=-2*r[i]*(t*dfdi[2+2*i]+imag(a)*dfdi[2+2*i+1]/w); + } + int i=0; + for(int k=0;k<3;k++) + *dfdx++=op(*dfdx,dfn[i++]); + for(int k=0;k<3;k++) + *dfdp++=op(*dfdp,dfn[i++]); + *dfdw=op(*dfdw,dfn[i++]); + *dfdpw=op(*dfdpw,dfn[i++]); + } + + + ///\en Compares the wave packet to another on a per component basis. + /// \return \retval 0 if all component differences are 0 within tolerance \a tol (EQUAL), + /// \retval -1 for LESS + /// \retval 2 for GREATER + int compare(const WavePacket &other, double tol=0.) const { + int res=compare_compl(a,other.a, tol); + if(res) + return res; + for(int i=0;i<3;i++){ + res=compare_compl(b[i],other.b[i]); + if(res) + return res; + } + return 0; + } + +}; + + /*double w=wk.get_width(); + Vector_3 r=wk.get_r(); + double t=3/(2*w*w*w); + fe_w[ic1+k1]+= t*E_der[s1][indw1+8*k1]+imag(wk.a)*E_der[s1][indw1+8*k1+1]/w; + fe_pw[ic1+k1]+=E_der[s1][indw1+8*k1+1]/(2*w*h_plank); + for(int i=0;i<3;i++){ + fe_x[ic1+k1][i]+= -2*real(wk.a)*E_der[s1][indw1+8*k1+2+2*i]-2*imag(wk.a)*E_der[s1][indw1+8*k1+2+2*i+1]; + fe_p[ic1+k1][i]+= (-E_der[s1][indw1+8*k1+2+2*i+1])*(m_electron/h_plank); //*(h_plank/m_electron); + fe_pw[ic1+k1]+=(r[i]*E_der[s1][indw1+8*k1+2+2*i+1]/w)/h_plank; + fe_w[ic1+k1]+=2*r[i]*(t*E_der[s1][indw1+8*k1+2+2*i]+imag(wk.a)*E_der[s1][indw1+8*k1+2+2*i+1]/w); + }*/ + + +/// constructs a conjugate packet +inline WavePacket conj(const WavePacket &wp){ + return WavePacket(conj(wp.a),conj(wp.b),conj(wp.lz)); +} + + +# endif // WAVEPACKET_H \ No newline at end of file diff --git a/lib/awpmd/ivutils/src/logexc.cpp b/lib/awpmd/ivutils/src/logexc.cpp new file mode 100644 index 0000000000..d29745fda1 --- /dev/null +++ b/lib/awpmd/ivutils/src/logexc.cpp @@ -0,0 +1,25 @@ +//# ifdef USE_STDAFX +//# include "stdafx.h" +//# endif + +# include "logexc.h" + +message_logger std_log; + +message_logger &message_logger::global(){ + if(!glogp){ + std_log.set_global(true); + } + return *glogp; +} + +message_logger *message_logger::glogp=NULL; +stdfile_logger default_log("",0,stdout,stderr,vblALLBAD|vblMESS1,vblFATAL,1); + +const char *fmt(const char *format,...){ + va_list args; + va_start(args,format); + static char buff[1024]; + vsnprintf(buff,1024,format,args); + return buff; +} diff --git a/lib/awpmd/license/COPYING b/lib/awpmd/license/COPYING new file mode 100644 index 0000000000..623b6258a1 --- /dev/null +++ b/lib/awpmd/license/COPYING @@ -0,0 +1,340 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. diff --git a/lib/awpmd/license/awpmd_license.txt b/lib/awpmd/license/awpmd_license.txt new file mode 100644 index 0000000000..346b4efbd8 --- /dev/null +++ b/lib/awpmd/license/awpmd_license.txt @@ -0,0 +1,8 @@ +AWPMD library license, Version 1 + +Copyright (c) 1992-2011 Ilya Valuev + +Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed. + +AWPMD library is free software, you can redistribute and modify it under the +terms of WXWINDOWS LIBRARY LICENSE, Version 3 (www.wxwidgets.org) diff --git a/lib/awpmd/license/wx_license.txt b/lib/awpmd/license/wx_license.txt new file mode 100644 index 0000000000..c91deed0bc --- /dev/null +++ b/lib/awpmd/license/wx_license.txt @@ -0,0 +1,53 @@ + wxWindows Library Licence, Version 3 + ==================================== + + Copyright (c) 1998 Julian Smart, Robert Roebling et al + + Everyone is permitted to copy and distribute verbatim copies + of this licence document, but changing it is not allowed. + + WXWINDOWS LIBRARY LICENCE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + This library is free software; you can redistribute it and/or modify it + under the terms of the GNU Library General Public Licence as published by + the Free Software Foundation; either version 2 of the Licence, or (at + your option) any later version. + + This library is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library + General Public Licence for more details. + + You should have received a copy of the GNU Library General Public Licence + along with this software, usually in a file named COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, + Boston, MA 02111-1307 USA. + + EXCEPTION NOTICE + + 1. As a special exception, the copyright holders of this library give + permission for additional uses of the text contained in this release of + the library as licenced under the wxWindows Library Licence, applying + either version 3 of the Licence, or (at your option) any later version of + the Licence as published by the copyright holders of version 3 of the + Licence document. + + 2. The exception is that you may use, copy, link, modify and distribute + under the user's own terms, binary object code versions of works based + on the Library. + + 3. If you copy code from files distributed under the terms of the GNU + General Public Licence or the GNU Library General Public Licence into a + copy of this library, as this licence permits, the exception does not + apply to the code that you add in this way. To avoid misleading anyone as + to the status of such modified files, you must delete this exception + notice from such code and/or adjust the licensing conditions notice + accordingly. + + 4. If you write modifications of your own for this library, it is your + choice whether to permit this exception to apply to your modifications. + If you do not wish that, you must delete the exception notice from such + code and/or adjust the licensing conditions notice accordingly. + + diff --git a/lib/awpmd/systems/interact/TCP/tcpdefs.h b/lib/awpmd/systems/interact/TCP/tcpdefs.h new file mode 100644 index 0000000000..e1e0b3553c --- /dev/null +++ b/lib/awpmd/systems/interact/TCP/tcpdefs.h @@ -0,0 +1,17 @@ +# ifndef _TCPDEFS_H +# define _TCPDEFS_H + +//Constants +const double sqr_2_over_3 = 0.816496580927726; +const double two_over_sqr_pi = 1.12837916709551; + +# ifndef PLASMA_UNITS // using GridMD units: A, eV, m_proton=1 +const double h_plank=0.06442021524615668; +const double h_sq=h_plank*h_plank; +const double m_electron=1./1836.1527556560675; +const double m_proton=1.; +const double coul_pref=14.39965172693122; // e^2/(4*pi*eps0), eV*A +const double eV_to_K=11604.447517053462; // Temperature in K correspondent to 1eV +# endif + +# endif \ No newline at end of file diff --git a/lib/awpmd/systems/interact/TCP/wpmd.cpp b/lib/awpmd/systems/interact/TCP/wpmd.cpp new file mode 100644 index 0000000000..ad5d8e26bc --- /dev/null +++ b/lib/awpmd/systems/interact/TCP/wpmd.cpp @@ -0,0 +1,890 @@ +# include "wpmd.h" + + + +// Calculates derivative overlap matrix IDD +void OverlapDeriv::calc_der_overlap(bool self, cdouble cc1, cdouble c2){ + cVector_3 I3 = I1 * ((bb_4a + 2.5) / w12.a); + cdouble I4 = I0 * ( bb_4a *(bb_4a + 5.) + 3.75 ) / w12.a / w12.a; + + // calculate derivatives <(phi_k)'_q_k | (phi_l)'_q_l>: + IDD.set(0, 0, I4 - (d1.l + d2.l)*I2 + d1.l*d2.l*I0 ); // over a_k_re and a_l_re + IDD.set(0, 1, i_unit*( I4 - (d1.l + d2.m)*I2 + d1.l*d2.m*I0 ) ); // over a_k_re and a_l_im + if(!self) + IDD.set(1, 0, i_unit1*( I4 + (d1.m - d2.l)*I2 - d1.m*d2.l*I0 ) ); // over a_k_im and a_l_re + else + IDD.set(1,0, conj(IDD(0,1))); + IDD.set(1, 1, I4 + (d1.m - d2.m)*I2 - d1.m*d2.m*I0 ); // over a_k_im and a_l_im + + for(int i=0;i<3;i++){ + IDD.set(0, (i+1)*2, -I3[i] + d1.l*I1[i] + d2.u[i]*(d1.l*I0 - I2) ); // over a_k_re and b_l_re + IDD.set(0, (i+1)*2+1, i_unit1*( I3[i] - d1.l*I1[i] + d2.v[i]*(I2 - d1.l*I0) ) ); // over a_k_re and b_l_im + IDD.set(1, (i+1)*2, i_unit *( I3[i] + d1.m*I1[i] + d2.u[i]*(I2 + d1.m*I0) ) ); // over a_k_im and b_l_re + IDD.set(1, (i+1)*2+1, -I3[i] - d1.m*I1[i] - d2.v[i]*(d1.m*I0 + I2) ); // over a_k_im and b_l_im + if(!self) { + IDD.set((i+1)*2, 0, -I3[i] + d2.l*I1[i] + d1.u[i]*(d2.l*I0 - I2) ); // over b_k_re and a_l_re + IDD.set((i+1)*2+1, 0, i_unit *( I3[i] - d2.l*I1[i] - d1.v[i]*(I2 - d2.l*I0) ) ); // over b_k_im and a_l_re + IDD.set((i+1)*2, 1, i_unit1*( I3[i] - d2.m*I1[i] + d1.u[i]*(I2 - d2.m*I0) ) ); // over b_k_re and a_l_im + IDD.set((i+1)*2+1, 1, -I3[i] + d2.m*I1[i] - d1.v[i]*(d2.m*I0 - I2) ); // over b_k_im and a_l_im + } + else{ + IDD.set((i+1)*2, 0, conj(IDD(0,(i+1)*2)) ); // over b_k_re and a_l_re + IDD.set((i+1)*2+1, 0, conj(IDD(0,(i+1)*2+1)) ); // over b_k_im and a_l_re + IDD.set((i+1)*2, 1, conj(IDD(1,(i+1)*2)) ); // over b_k_re and a_l_im + IDD.set((i+1)*2+1, 1, conj(IDD(1,(i+1)*2+1)) ); // over b_k_im and a_l_im + } + + for(int j=0;j<3;j++){ + if(!self || j>=i){ + cdouble I2ij = I0 / w12.a * + (i==j ? w12.b[i]*w12.b[i] / w12.a / 4 + 0.5 + : w12.b[i]*w12.b[j] / w12.a / 4); + // over b_k_re and b_l_re + IDD.set((j+1)*2, (i+1)*2, I2ij + d1.u[i]*I1[j] + d2.u[j]*(I1[i] + d1.u[i]*I0) ); + // over b_k_re and b_l_im + IDD.set((j+1)*2, (i+1)*2+1, i_unit *( I2ij + d1.u[i]*I1[j] + d2.v[j]*(I1[i] + d1.u[i]*I0) ) ); + // over b_k_im and b_l_re + if(!self || i!=j) + IDD.set((j+1)*2+1, (i+1)*2, i_unit1*( I2ij - d1.v[i]*I1[j] + d2.u[j]*(I1[i] - d1.v[i]*I0) ) ); + else + IDD.set((j+1)*2+1, (i+1)*2, conj(IDD((i+1)*2,(j+1)*2+1))); + // over b_k_im and b_l_im + IDD.set((j+1)*2+1,(i+1)*2+1, I2ij - d1.v[i]*I1[j] + d2.v[j]*(I1[i] - d1.v[i]*I0) ); + } + else{ // self && j0) + w=w0; + pw=0.; + } + + double rw; + if(Lextra>0){ // width PBC, keeping the width are within [0,Lextra] + w=fmod(w,Lextra); + if(w<0) w+=Lextra; + rw=w; // WP width for energy evaluation is within [0, L/2] + if(rw > Lextra/2) rw = Lextra - rw; + } + else + rw=w; + + WavePacket wp; + wp.init(rw,x,v*mass*one_h,pw*one_h); + return wp; +} + + +void AWPMD::resize(int flag){ + for(int s=0;s<2;s++){ + //0. resizing matrices + Y[s].init(ne[s],1); + O[s].init(ne[s],1); + Oflg[s].init(ne[s],1); + //Te[s].init(nel,1); + //Tei[s].init(nel,1); + Eep[s].assign((size_t)nwp[s],0); + Eeip[s].assign((size_t)nwp[s],0); + Eeep[s].assign((size_t)nwp[s],0); + Ewp[s].assign((size_t)nwp[s],0); + + if(flag&(0x8|0x4) && approx!=HARTREE){ //electron forces, L and M are needed + M[s].init(ne[s],nvar[s]); + L[s].init(ne[s],nvar[s]); + } + } + Eiep.assign((size_t)ni,0); + Eiip.assign((size_t)ni,0); + + +} + + +//e sets Periodic Boundary Conditions +//e using bit flags: 0x1 -- PBC along X +//e 0x2 -- PBC along Y +//e 0x4 -- PBC along Z +//e cell specifies the lengths of the simulation box in all directions +//e if PBCs are used, the corresponding coordinates of electrons and ions +//e in periodic directions must be within a range [0, cell[per_dir]) +//e @returns 1 if OK +int AWPMD::set_pbc(const Vector_3P pcell, int pbc_){ + if(!pcell) + pbc=0; + else{ + pbc=pbc_; + cell=*pcell; + } + return 1; +} + +//e setup elctrons: forms internal wave packet representations +//e if PBCs are used the coords must be within a range [0, cell) +int AWPMD::set_electrons(int s, int n, Vector_3P x, Vector_3P v, double* w, double* pw, double mass, double *q) +{ + if(s < 0 || s > 1) + return LOGERR(-1,fmt("AWPMD.set_electrons: invaid s setting (%d)!",s),LINFO); + + norm_matrix_state[s] = NORM_UNDEFINED; + nwp[s]=ne[s]=n; + nvar[s]=8*n; + wp[s].resize(n); + + partition1[s].clear(); + for(int i=0;i0){ // width PBC + if(w1>Lextra-w1){ + w1=-(Lextra-w1); // '-' is to change derivative sign + sgn1=-1; + } + }*/ + Vector_3 r1=wp[s1][c1].get_r(); + Vector_3 p=wp[s1][c1].get_p()*h_plank; + Vector_3 pw=wp[s1][c1].get_pwidth()*h_plank; + // energy contribution + Ee[s1] += (p.norm2()+pw.norm2())/(2*me); + Ew += h2_me*9./(8.*w1*w1); + if(constraint == HARM) Ew += harm_w0_4 * w1*w1; + // width force contribution + //double dE=2*Epot/w; + //if(d->fw1)d->fw1[c1]+=dE; + //if(fw2 && fw2!=fw1)fw2[c1]+=dE; + + // e-e interaction + for(int s2=s1;s2<2;s2++){ + for(int c2=(s1==s2 ? c1+1 : 0) ;c20){ // width PBC + if(w2>Lextra-w2){ + w2=-(Lextra-w2); // '-' is to change derivative sign + sgn2=-1; + } + }*/ + double wsq=w1*w1+w2*w2; + double argw=sqrt((2./3.)*wsq); + + //double arg=r12/argw; + //double erfa=erf(arg); + double Epot=coul_pref*erf_div(r12,1./argw); //erfa/r12; + Eee+=Epot; + + // force contribution + /*double dEw=coul_pref*two_over_sqr_pi*exp(-arg*arg)/argw; + double dEpot=(Epot-dEw)/r12; + if(!d->fixw){ + dEw/=wsq; + if(d->fw1 && c1>=0){ + d->fw1[c1]+=sgn1*dEw*w1; + } + if(d->fw2){ + d->fw2[c2]+=sgn2*dEw*w2; + } + }*/ + } + } + // e-i interaction + double wsq1=w1*w1; + double argw=sqr_2_over_3*w1; + for(int i=0;ifixw){ + dEw/=wsq; + if(d->fw1 && c1>=0){ + d->fw1[c1]+=sgn1*dEw*w1; + } + }*/ + } + } + } + if(calc_ii) + interaction_ii(flag,fi); + return 1; +} + + +//e initializes internal buffers for calculations (set_electrons must be called first) +//int init(){} + +//e calculates interaction in the system of ni ions + electrons +//e the electonic subsystem must be previously setup by set_electrons, ionic by set_ions +//e the iterators are describing ionic system only +// 0x1 -- give back ion forces +// 0x2 -- add ion forces to the existing set +// 0x4 -- calculate derivatives for electronic time step (NOT IMPLEMENTED) +//e if PBCs are used the coords must be within a range [0, cell) +int AWPMD::interaction(int flag, Vector_3P fi, Vector_3P fe_x, + Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){ + if(approx==HARTREE) + return interaction_hartree(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c); + // 0. resizing the arrays if needed + resize(flag); + // 1. clearing forces + clear_forces(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c); + + //2. calculating overlap matrix + for(int s=0;s<2;s++){ + int nes = ne[s]; + if(nes == 0) continue; + + for(int k=0;k ovl_tolerance); + } + } + O[s] = Y[s]; // save overlap matrix + + //3. inverting the overlap matrix + int info=0; + if(nes){ + /*FILE *f1=fopen(fmt("matrO_%d.d",s),"wt"); + fileout(f1,Y[s],"%15g"); + fclose(f1);8*/ + + ZPPTRF("L",&nes,Y[s].arr,&info); + // analyze return code here + if(info<0) + return LOGERR(info,fmt("AWPMD.interacton: call to ZPTRF failed (exitcode %d)!",info),LINFO); + ZPPTRI("L",&nes,Y[s].arr,&info); + if(info<0) + return LOGERR(info,fmt("AWPMD.interacton: call to ZPTRI failed (exitcode %d)!",info),LINFO); + + + /*f1=fopen(fmt("matrY_%d.d",s),"wt"); + fileout(f1,Y[s],"%15g"); + fclose(f1);*/ + } + + // Clearing matrices for electronic forces + if(flag&0x4){ + Te[s].Set(0.); + Tei[s].Set(0.); + } + } + + Vector_3 ndr; + // calculating single particle contribution + for(int s=0;s<2;s++){ + Ee[s]=Eei[s]=0.; + for(int k=0;k1e-10){ + cdouble arg=ngkli*c; + cdouble dEw=-coul_pref*qi[i]*I0kl*two_over_sqr_pi*exp(-arg*arg)*c; + dE=(dE-dEw)/ngkli; + dE*=Y[s](l,k); + Vector_3 dir=-real(gkli); + dir.normalize(); + fi[i]+=(l==k ? 1. : 2.)*real(dE)*dir; + } + } + } + dE=sum; + if(flag&0x4) // matrix needed only for electronic forces + Tei[s].set(k,l,dE); + // energy component (trace) + dE*=Y[s](l,k); + Eei[s]+=(l==k ? 1. : 2.)*real(dE); + } + } + } + + // calculating e-e interaction + Eee = Ew = 0.; + // same spin + for(int s=0;s<2;s++){ // spin + for(int k=0;k& Norms = Norm[s]; + chmatrix& Ys = Y[s]; + smatrix& Oflgs = Oflg[s]; + + // Allocate of vectors and matrices + Norms.init(nes8,1); + IDD.init(nes8,1); + if(ID.size() != nnes8) + ID.resize(nnes8), IDYs.resize(nnes8), ipiv.resize(nes8); + + // Calculate first and second derivatives + for(k=0;k -mu, v -> -v !!! + + for(l=0;l: + int idx = k + l*nes8; + if(k != l) { + ID[idx] = dl.l*I0 - I2; // over a_l_re + ID[idx+nes] = i_unit*(dl.m*I0 - I2); // over a_l_im + for(i=0;i<3;i++){ + ID[idx+((i+1)*2)*nes] = dl.u[i]*I0 + I1[i]; // over b_l_re + ID[idx+((i+1)*2+1)*nes] = i_unit*(dl.v[i]*I0 + I1[i]); // over b_l_im + } + } else { // k == l + ID[idx] = i_unit*imag(dl.l); // over a_l_re + ID[idx+nes] = i_unit*(dl.m - I2); // over a_l_im + for(i=0;i<3;i++){ + ID[idx+((i+1)*2)*nes] = dl.u[i] + I1[i]; // over b_l_re + ID[idx+((i+1)*2+1)*nes] = 0.; // over b_l_im + } + } + + if(k <= l) { + cVector_3 I3 = I1 * ((bb_4a + 2.5) / wkl.a); + cdouble I4 = I0 * ( bb_4a *(bb_4a + 5.) + 3.75 ) / wkl.a / wkl.a; + + // calculate derivatives <(phi_k)'_q_k | (phi_l)'_q_l>: + IDD.set(k8, l8, I4 - (dk.l + dl.l)*I2 + dk.l*dl.l*I0 ); // over a_k_re and a_l_re + IDD.set(k8, l8+1, i_unit*( I4 - (dk.l + dl.m)*I2 + dk.l*dl.m*I0 ) ); // over a_k_re and a_l_im + if(k != l) IDD.set(k8+1, l8, i_unit1*( I4 + (dk.m - dl.l)*I2 - dk.m*dl.l*I0 ) ); // over a_k_im and a_l_re + IDD.set(k8+1, l8+1, I4 + (dk.m - dl.m)*I2 - dk.m*dl.m*I0 ); // over a_k_im and a_l_im + + for(i=0;i<3;i++){ + IDD.set(k8, l8+(i+1)*2, -I3[i] + dk.l*I1[i] + dl.u[i]*(dk.l*I0 - I2) ); // over a_k_re and b_l_re + IDD.set(k8, l8+(i+1)*2+1, i_unit1*( I3[i] - dk.l*I1[i] + dl.v[i]*(I2 - dk.l*I0) ) ); // over a_k_re and b_l_im + IDD.set(k8+1, l8+(i+1)*2, i_unit *( I3[i] + dk.m*I1[i] + dl.u[i]*(I2 + dk.m*I0) ) ); // over a_k_im and b_l_re + IDD.set(k8+1, l8+(i+1)*2+1, -I3[i] - dk.m*I1[i] - dl.v[i]*(dk.m*I0 + I2) ); // over a_k_im and b_l_im + if(k != l) { + IDD.set(k8+(i+1)*2, l8, -I3[i] + dl.l*I1[i] + dk.u[i]*(dl.l*I0 - I2) ); // over b_k_re and a_l_re + IDD.set(k8+(i+1)*2+1, l8, i_unit *( I3[i] - dl.l*I1[i] - dk.v[i]*(I2 - dl.l*I0) ) ); // over b_k_im and a_l_re + IDD.set(k8+(i+1)*2, l8+1, i_unit1*( I3[i] - dl.m*I1[i] + dk.u[i]*(I2 - dl.m*I0) ) ); // over b_k_re and a_l_im + IDD.set(k8+(i+1)*2+1, l8+1, -I3[i] + dl.m*I1[i] - dk.v[i]*(dl.m*I0 - I2) ); // over b_k_im and a_l_im + } + + for(j=0;j<3;j++){ + cdouble I2ij = I0 / wkl.a * + (i==j ? wkl.b[i]*wkl.b[i] / wkl.a / 4 + 0.5 + : wkl.b[i]*wkl.b[j] / wkl.a / 4); + // over b_k_re and b_l_re + IDD.set(k8+(j+1)*2, l8+(i+1)*2, I2ij + dk.u[i]*I1[j] + dl.u[j]*(I1[i] + dk.u[i]*I0) ); + // over b_k_re and b_l_im + IDD.set(k8+(j+1)*2, l8+(i+1)*2+1, i_unit *( I2ij + dk.u[i]*I1[j] + dl.v[j]*(I1[i] + dk.u[i]*I0) ) ); + // over b_k_im and b_l_re + if(k != l) IDD.set(k8+(j+1)*2+1, l8+(i+1)*2, i_unit1*( I2ij - dk.v[i]*I1[j] + dl.u[j]*(I1[i] - dk.v[i]*I0) ) ); + // over b_k_im and b_l_im + IDD.set(k8+(j+1)*2+1, l8+(i+1)*2+1, I2ij - dk.v[i]*I1[j] + dl.v[j]*(I1[i] - dk.v[i]*I0) ); + } // j + } // i + } // if(k <= l) + } // k + } // l + + // Calculate matrix product IDYs_(k,q_j) = Ys_(k,l) * + for(qj=0; qj * IDYs_(k,q_j) + cdouble sum = 0.; + for(k=0;k::iterator mi=Norms.fix_first(8*i+k,0); + for(j=i ;j(mi+8*j,mi+8*j,mi+8*j+3,mi+8*j+6,mi+8*j+7); + } + }// finished line of blocks from right + for(k= 8*i;k::iterator mi=Norms.fix_second(8*i,k); + wi.int2phys_der< eq_second >(mi,mi,mi+3,mi+6,mi+7); + }// finished line of blocks from left + + for(k=0;k<8;k++){ // filling the lower triangle according to antisymmetry + for(j=8*i+8;j::iterator mi=Norms.fix_first(8*i+k,8*j); + wj.int2phys_der< eq_second >(mi,mi,mi+3,mi+6,mi+7); + } + for(k=0;k<8;k++){ + // iterator to list all N(8*i+k,*) with fixed 8*i+k + sqmatrix::iterator mi=Norms.fix_second(8*i,8*j+k); + wi.int2phys_der< eq_second >(mi,mi,mi+3,mi+6,mi+7); + } + if(i!=j){ + for(int k1=0;k1<8;k1++){ + for(int k2=0;k2<8;k2++) + Norms(8*j+k1,8*i+k2)=-Norms(8*i+k2,8*j+k1); + } + } + } + } +# endif + + norm_matrix_state[s] = NORM_CALCULATED; +} + +//e Norm matrix LU-factorization +void AWPMD::norm_factorize(int s) { + if( norm_matrix_state[s] != NORM_CALCULATED) norm_matrix(s); + + int nes8 = ne[s]*8, info; + DGETRF(&nes8, &nes8, Norm[s].arr, &nes8, &ipiv[0], &info); + if(info < 0) + LOGERR(info,fmt("AWPMD.norm_factorize: call to DGETRF failed (exitcode %d)!",info),LINFO); + + norm_matrix_state[s] = NORM_FACTORIZED; +} + + +//e Norm matrix inversion +void AWPMD::norm_invert(int s) { + if( norm_matrix_state[s] != NORM_FACTORIZED) norm_factorize(s); + + int nes8 = ne[s]*8, info; + int IDD_size = (int)IDD.get_datasize(nes8); + + DGETRI(&nes8, Norm[s].arr, &nes8, &ipiv[0], (double*)IDD.arr, &IDD_size, &info); // use IDD for work storage + if(info < 0) + LOGERR(info,fmt("AWPMD.norm_invert: call to DGETRI failed (exitcode %d)!",info),LINFO); + + norm_matrix_state[s] = NORM_INVERTED; +} + + +//e Get the determinant of the norm-matrix for the particles with spin s +double AWPMD::norm_matrix_det(int s) { + double det = 1.; + int nes8 = ne[s]*8; + + if(!nes8) return det; + if(norm_matrix_state[s] != NORM_FACTORIZED) norm_factorize(s); + + sqmatrix& Norms = Norm[s]; + for(int i=0; i& Norms = Norm[s]; + for(int i=0; i1) + return -1; // invalid spin: return LOGERR(-1,fmt("AWPMD.get_electrons: invaid spin setting (%d)!",spin),LINFO); + if(mass<0) + mass=me; + for(int i=0;i, some fixes to UHF + * + * Revision 1.35 2011/05/25 21:30:48 morozov + * Compatibility with ICC 11.1 + * + * Revision 1.34 2011/05/25 05:23:43 valuev + * fixed variable transformation for norm matrix + * + * Revision 1.33 2011/05/24 19:54:32 valuev + * fixed sqmatrix::iterator + * + * Revision 1.32 2011/05/21 23:06:49 valuev + * Norm matrix transform to pysical variables + * + * Revision 1.31 2011/05/20 21:39:49 valuev + * separated norm calculation + * + * Revision 1.30 2011/05/14 18:56:19 valuev + * derivative for ee split interactions + * + * Revision 1.29 2011/05/05 08:56:02 valuev + * working split wp version + * + * Revision 1.28 2011/05/04 16:48:52 valuev + * fixed syntax + * + * Revision 1.27 2011/05/04 09:04:48 valuev + * completed wp_split (except for ee forces) + * + * Revision 1.26 2011/04/29 03:07:20 valuev + * new split wp features + * + * Revision 1.25 2011/04/22 09:52:49 valuev + * working on split WP + * + * Revision 1.24 2011/04/20 08:43:09 valuev + * started adding split packet WPMD + * + * Revision 1.23 2010/09/03 12:17:48 morozov + * The order of parameters in Norm matrix is changed to mimic the order of the single WP parameter storage + * + * Revision 1.22 2009/08/27 00:01:36 morozov + * First working MPI equilibration + * + * Revision 1.21 2009/04/14 14:44:10 valuev + * fixed momentum calculation in hartree model, added "fix" constraint and model="hartree" + * to parameters + * + * Revision 1.20 2009/04/13 17:00:45 morozov + * Fixed norm-matrix ratio in AWPMC algorithm + * + * Revision 1.19 2009/04/06 17:00:28 morozov + * Fixed Hartree version of WPMC + * + * Revision 1.18 2009/04/01 10:06:37 valuev + * added Hartee factorization to AWPMD + * + * Revision 1.17 2009/03/24 16:10:05 morozov + * Fixed errors in Norm-matrix calculation related to PBC + * + * Revision 1.16 2009/03/17 11:40:04 morozov + * The prefactor of NormMatrix is corrected + * + * Revision 1.15 2008/07/23 16:42:12 valuev + * Added AWPMD Monte-Carlo + * + * Revision 1.14 2008/07/23 15:58:32 valuev + * *** empty log message *** + * + * Revision 1.13 2008/07/21 02:23:22 morozov + * *** empty log message *** + * + * Revision 1.12 2008/07/18 18:15:31 morozov + * *** empty log message *** + * + * Revision 1.11 2008/05/29 13:33:05 valuev + * VASP band structure + * + * Revision 1.10 2008/05/14 17:17:26 morozov + * Passed 2- and 3-electron test. Added Norm matrix. + * + * Revision 1.9 2008/05/05 17:29:32 morozov + * Fixed errors with Hermitian matrix indeces. Redesigned cVector_3 operations. + * + * Revision 1.8 2008/04/28 22:16:45 valuev + * restructured coulomb term + * + * Revision 1.7 2008/04/28 09:54:13 valuev + * corrected summation for Eee part + * +*******************************************************************************/ +# ifndef WPMD_H +# define WPMD_H + +/** @file wpmd.h + @brief Classes for Wave Packet Molecular Dynamics of two component plasma. */ + +# ifndef _USE_MATH_DEFINES +# define _USE_MATH_DEFINES +# endif +# include +# include +# include +# include "logexc.h" +# include "cvector_3.h" +# include "pairhash.h" +# include "TCP/tcpdefs.h" +# include "wavepacket.h" +# include "erf.h" +# include "cerf.h" + + +using namespace std; + +# include "lapack_inter.h" + +typedef hmatrix chmatrix; + +const cdouble i_unit = cdouble(0.,1.), i_unit1 = -i_unit; +const double h_plank2 = h_plank * 2.; + +//cdouble ccerf(const cdouble &a); + +//e calculates cerf(c*z)/z +inline cdouble cerf_div(const cdouble &z, const cdouble &c=i_unit){ + if((fabs(real(z))+fabs(imag(z)))<1e-8) + return c*two_over_sqr_pi; + else + return cerf(z*c)/z; +} + +//e calculates erf(c*z)/z +inline double erf_div(const double &z, double c=1){ + if(fabs(z)<1e-8) + return c*two_over_sqr_pi; + else + return erf(z*c)/z; +} + +template +struct _myrefpair{ + T1& first; + T2& second; + _myrefpair(T1& a, T2 &b):first(a),second(b){} + _myrefpair operator=(const pair &other){ + first=other.first; + second=other.second; + return *this; + } +}; + +template< class T1, class T2> +_myrefpair _mytie(T1& var1, T2 &var2){ + return _myrefpair(var1, var2); +} + +inline pair operator*(const pair &right, double left){ + return make_pair(right.first*left,right.second*left); +} + +// Auxilary class to handle the normalizing term derivatives +class NormDeriv +{ +public: + cdouble l; // lambda = (f over a_re) + double m; // mu = (f over a_im) / i + cVector_3 u; // u = (f over b_re) + Vector_3 v; // v = (f over b_im) / i + + NormDeriv() {} + NormDeriv(const WavePacket& wp) { set(wp); } + + //e Create NormDeriv object and calculate the derivatived for the given WP + void set(const WavePacket& wp){ + Vector_3 br = real(wp.b), bi = imag(wp.b); + double ar = real(wp.a), ai = imag(wp.a); + double i_2ar = 0.5 / ar, ai_ar = ai / ar; + + v = (-i_2ar) * br; + m = v.norm2(); + u = v * (i_unit1 * ai_ar) - wp.b * i_2ar; + l = (1.5*i_2ar + m) + cdouble(0.,2.) * ( (br*bi)*i_2ar*i_2ar - m*ai_ar ); + } +}; + +inline NormDeriv conj(const NormDeriv& src){ + NormDeriv dst; + dst.l = conj(src.l); + dst.m = -src.m; + dst.u = conj(src.u); + dst.v = - src.v; + return dst; +} + +///\en Auxilary class to handle derivatives of overlaps +class OverlapDeriv{ +public: + WavePacket w1, w2, w12; + NormDeriv d1, d2; + cdouble I0, I2; + cVector_3 I1; + cdouble bb_4a; + sqmatrix IDD; + + + OverlapDeriv():I0(0),I1(0),IDD(10){} + + void set1(const WavePacket& w1_) { + w1=w1_; + d1.set(w1); + d1=conj(d1); + } + + //e Create NormDeriv object and calculate the derivatived for the given WP + void set2(const WavePacket& w2_, const cdouble *I0_=NULL){ + w2=w2_; + d2.set(w2); + w12=conj(w1)*w2; + if(!I0_) + I0 = w12.integral(); + else + I0=*I0_; + I1 = w12.b * (I0 / w12.a / 2); + bb_4a = w12.b.norm2() / w12.a / 4; + I2 = I0 * (bb_4a + 1.5) / w12.a; + } + + cdouble da2_re() const { + return (d2.l*I0 - I2); + } + + cdouble da2_im() const { + return i_unit*(d2.m*I0 - I2); + } + + cdouble da1_re() const { + return (d1.l*I0 - I2); + } + + cdouble da1_im() const { + return -i_unit*(d1.m*I0 - I2); + } + + cdouble db2_re(int i) const { + return d2.u[i]*I0 + I1[i]; + } + + cdouble db2_im(int i) const { + return i_unit*(d2.v[i]*I0 + I1[i]); + } + + cdouble db1_re(int i) const { + return d1.u[i]*I0 + I1[i]; + } + + cdouble db1_im(int i) const { + return -i_unit*(d1.v[i]*I0 + I1[i]); + } + + ///\en Calculates derivative overlap matrix IDD + void calc_der_overlap(bool self=false, cdouble cc1=0., cdouble c2=0.); + + +}; + + +class AWPMD { +//protected: +public: + int ne[2], ni; + int nwp[2]; ///<\en number of wavepackets (same as ne for unsplit version) + int nvar[2]; ///<\en full number of dynamic variables for each spin + chmatrix O[2], Y[2], Te[2], Tei[2]; + smatrix Oflg[2]; ///<\en equals 0 for non-overlaping packets + sqmatrix Norm[2]; ///<\en Norm matrix + vector wp[2]; ///<\en wave packets for electrons (spin up =0 and down =1) + vector qe[2]; ///<\en electron charges + vector qi; ///<\en ion charges + vector xi; ///<\en ion coordinates + int pbc; ///<\en pbc flag + Vector_3 cell; ///<\en cell coordinates (L1,L2,L3) + double Lextra; ///<\en width PBC, unset if negative + double harm_w0_4; + double w0; + int calc_ii; ///<\en flag indicating whether to calculate ion-ion interaction + int norm_needed; ///<\en flag indicating whether to prepare norm matrix data in interaction function + + enum {NORM_UNDEFINED, NORM_CALCULATED, NORM_FACTORIZED, NORM_INVERTED}; + int norm_matrix_state[2]; + + // Arrays for temporal data + chmatrix IDD; // Second derivatives of the overlap integral (used in Norm matrix) + vector ID, IDYs; // First derivatives of the overlap integral (used in Norm matrix) + vector ipiv; // The pivot indices array + + recmatrix L[2]; ///<\en overlap derivative matrix for each spin + recmatrix M[2]; ///<\en 'normalized' overlap derivative matrix for each spin: M=YL + +public: + enum {NONE=0, HARM, FIX, RELAX} constraint; + + ///\em Sets approximation level for quantum problem: \n + /// HARTREE Hartree product (no antisymmetrization) \n + /// DPRODUCT product of det0*det1 of antisymmetrized functions for spins 0, 1 \n + /// UHF unrestricted Hartree-Fock + enum APPROX {HARTREE, DPRODUCT, UHF } approx; + ///\em Sets overlap matrix element to zero if the overlap norm is less than this value + double ovl_tolerance; + + double Ee[2]; ///<\en electron kinetic energy for each spin + double Eei[2]; ///<\en electron-ion energy for each spin + double Eee, Ew; ///<\en electron-electron energy + double Eii; ///<\en ion-ion energy + double Edk; ///<\en sum of diagonal kinetic energy terms + double Edc; ///<\en sum of diagonal Coulomb energy terms + + vector Eep[2]; ///<\en per particle electron kinetic energy for each spin + vector Eeip[2]; ///<\en per particle electron-ion energy for each spin + vector Eeep[2]; ///<\en per particle electron-electron energy for each spin + vector Ewp[2]; ///<\en per particle restrain energy for each spin + vector Eiep; ///<\en per particle ion-electron energy + vector Eiip; ///<\en per particle ion-ion energy + + + ///\en \{ Conversion constants that depend on the unit system used (for LAMMPS compatibility). + /// Default is GRIDMD units. Change them according to your unit system. + double me; ///<\en electron mass (LAMMPS: me in the appropriate unit system) + double one_h; ///<\en inverse of Plancks constant (LAMMPS: conversion [(m*v)/h] to [distance] ) + double h2_me; ///<\en Plancks constant squared divided by electron mass (LAMMPS: conversion [h^2/(m*r^2)] to [Energy] ) + double coul_pref; ///<\en Coulomb prefactor (e2 for GRIDMD) (LAMMPS: conversion [q^2/r] to [Energy] ) + /// \} + + ///\en 0 -- indicates that the inter-partition force should be full, and energy half,\n + /// 1 -- inter-partition force and energy counts one half (LAMMPS compatibility) + int newton_pair; + + //int myid; ///<\en id for partitions + + ///\en Partition arrays storing the tags of particles. The initial tags should be >0. + /// If the tag stored is <0, then the particle is ghost with -tag. + /// partition1[2] is for ions, 0, 1 for each electron spin + vector partition1[3]; + //vector partition2[3]; ///<\en 2 for ions + + + int tag_index(int i, int j){ + return i==j ? -1 : (i>j ? (i-2)*(i-1)/2+j : (j-2)*(j-1)/2+i ); + } + + + ///\en 1 -- all my, -1 all other, 2 -- my mixed term, -2 -- other mixed term + int check_ee(int s1,int icj1,int ick2){ + //printf(" (%d %d) ",partition1[s1][icj1],partition1[s1][ick2]); + int c1=(int)(partition1[s1][icj1]>0); + int c2=(int)(partition1[s1][ick2]>0); + int res; + if(c1!=c2){ // mixed + int tag1=abs(partition1[s1][icj1]); + int tag2=abs(partition1[s1][ick2]); + int num=tag_index(tag1-1,tag2-1); + if(num<0){ // compare wave packets + int cmp= s1<2 ? + wp[s1][icj1].compare(wp[s1][ick2],1e-15) : + compare_vec(xi[icj1],xi[ick2],1e-15); + if((cmp>0 && c1) || (cmp<0 && c2)) + res= 2; // my mixed term + else + res= -2; // not my term + } + else // parity check + res=num%2 ? 2 : -2; + } + else if(c1) + res=1; // all my + else + res=-1; // all other + return res; + } + + ///\en Returns electron-electron inter-partition multipliers for energy (first) and force (second) + /// for a 4- and 2- electron additive terms (all inter-partition interactions are + /// calculated only once based on particle tags) + /// If force multiplier is zero, then the term may be omitted (energy will also be zero). + /// NOW ASSIGNS BASED ON THE FIRST PAIR ONLY + pair check_part1(int s1,int icj1,int ick2, int s2=-1,int icj3=-1,int ick4=-1){ + int res=check_ee(s1,icj1,ick2); + if(res==1){ // my term + //printf(" *\n"); + return make_pair(1.,1.); // all at my partition + } + else if(res==-1){ + //printf(" \n"); + return make_pair(0.,0.); // all at other partition + } + else if(res==2){ + //printf(" *\n"); + return make_pair(1.,1.); // my inter-partition + } + else if(res==-2){ + //printf(" \n"); + return make_pair(0., newton_pair ? 0.0 : 1. ); // other inter-partition: must add force if newton comm is off + } + return make_pair(0.,0.); // nonsense + } + + ///\en Returns elctron-ion inter-partition multipliers for energy (first) and force (second) + /// for ion-electron additive terms (all inter-partition interactions are + /// calculated only once based on particle tags) + /// If force multiplier is zero, then the term may be omitted (energy will also be zero). + /// BASED ON ION ATTACHMENT + pair check_part1ei(int s1,int icj1,int ick2, int ion){ + //printf("%d ",partition1[2][ion]); + int ci=(int)(partition1[2][ion]>0); + + if(!newton_pair){ // care about mixed terms + int cee=check_ee(s1,icj1,ick2); + if((cee==2 || cee==-2) || (ci && cee==-1) || (!ci && cee==1)) // all mixed variants + make_pair(0., 1. ); // other inter-partition: must add force if newton comm is off + } + if(ci){ + //printf(" *\n"); + return make_pair(1.,1.); // my term + } + else{ + //printf(" \n"); + return make_pair(0.,0.); // all at other partition + } + } + + ///\en Returns ion-ion inter-partition multipliers for energy (first) and force (second) + /// for ion-electron additive terms (all inter-partition interactions are + /// calculated only once based on particle tags) + /// If force multiplier is zero, then the term may be omitted (energy will also be zero). + pair check_part1ii(int ion1, int ion2){ + return check_part1(2,ion1,ion2); + } + + + + AWPMD():pbc(0),Lextra(-1),constraint(NONE),newton_pair(1) { + nvar[0]=nvar[1]=ne[0]=ne[1]=ni=0; + norm_matrix_state[0] = norm_matrix_state[1] = NORM_UNDEFINED; + ovl_tolerance=0.; + approx = DPRODUCT; + + me=m_electron; + one_h=1./h_plank; + h2_me=h_sq/me; + coul_pref=::coul_pref; + + calc_ii=0; + norm_needed=0; + } + +protected: + //e translates wp2 to the nearest image postion relative to wp1 + //e gets the translation vector + Vector_3 move_to_image(const WavePacket &wp1, WavePacket &wp2) const { + Vector_3 r1=wp1.get_r(); + Vector_3 r2=wp2.get_r(); + Vector_3 dr=r2-r1; + Vector_3 ndr=dr.rcell(cell,pbc); // [0,L) + ndr=ndr.rcell1(cell,pbc); // [-L/2,L/2) + ndr-=dr; + wp2=wp2.translate(ndr); // wln.b=wln.b+ndr.a + return ndr; + } + + //e gets the overlap packet taking PBC into account + WavePacket pbc_mul(const WavePacket &wp1, const WavePacket &wp2) const { + if(!pbc) + return wp1*conj(wp2); + Vector_3 r1=wp1.get_r(); + Vector_3 r2=wp2.get_r(); + Vector_3 dr=r2-r1; // distance + Vector_3 drn=dr.rcell1(cell,pbc); // distance within PBC + Vector_3 rtrans=drn-dr; // new location of wp2 according to PBC (nearest image) + WavePacket wpn=wp2.translate(rtrans); + wpn=wp1*(conj(wpn)); + // reducing the result to elementary cell + //r1=wpn.get_r(); + //r2=r1.rcell(cell,pbc); + //dr=r2-r1; + //wpn=wpn.translate(dr); + return wpn; + } + ///\en resizes all internal arrays according to new electrons added + virtual void resize(int flag); +public: + + + + ///\en Prepares to setup a new system of particles using \ref add_ion() and add_electron(). + /// There is no need to call this function when using + /// \ref set_electrons() and \ref set_ions() to setup particles. + virtual void reset(){ + for(int s=0;s<2;s++){ + nwp[s]=ne[s]=nvar[s]=0; + wp[s].clear(); + qe[s].clear(); + partition1[s].clear(); + //partition2[s].clear(); + } + partition1[2].clear(); + ni=0; + xi.clear(); + qi.clear(); + } + + //e sets Periodic Boundary Conditions + //e using bit flags: 0x1 -- PBC along X + //e 0x2 -- PBC along Y + //e 0x4 -- PBC along Z + //e cell specifies the lengths of the simulation box in all directions + //e if PBCs are used, the corresponding coordinates of electrons and ions + //e in periodic directions must be within the range [0, cell[per_dir]) + //e @returns 1 if OK + int set_pbc(const Vector_3P pcell=NULL, int pbc_=0x7); + + ///\en Setup electrons: forms internal wave packet representations. + /// If PBCs are used the coords must be within a range [0, cell). + /// Default electron mass is AWPMD::me. + /// Default (q=NULL )electron charges are -1. + int set_electrons(int spin, int n, Vector_3P x, Vector_3P v, double* w, double* pw, double mass=-1, double *q=NULL); + + //e setup ion charges and coordinates + //e if PBCs are used the coords must be within a range [0, cell) + int set_ions(int n, double* q, Vector_3P x); + + ///\en Adds an ion with charge q and position x, + /// \return id of the ion starting from 0 + /// The tags must be nonzero, >0 for the local particle, <0 for ghost particle. + /// Unique particle id is abs(tag). + /// Default tag (0) means inserting the current particle id as local particle. + int add_ion(double q, const Vector_3 &x, int tag=0){ + qi.push_back(q); + xi.push_back(x); + ni=(int)xi.size(); + if(tag==0) + tag=ni; + partition1[2].push_back(tag); + return ni-1; + } + + + //e calculates interaction in the system of ni ions + electrons + //e the electonic subsystem must be previously setup by set_electrons, ionic by set_ions + //e the iterators are describing ionic system only + // 0x1 -- give back ion forces + // 0x2 -- add ion forces to the existing set + // 0x4 -- calculate derivatives for electronic time step (NOT IMPLEMENTED) + //e if PBCs are used the coords must be within a range [0, cell) + virtual int interaction(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL, + Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL); + + //e same as interaction, but using Hartee factorization (no antisymmetrization) + virtual int interaction_hartree(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL, + Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL); + + ///\en Calculates ion-ion interactions and updates Eii and ion forces if requested. This function + /// is called form intaraction() and interaction_hartree if calc_ii is set. + virtual int interaction_ii(int flag,Vector_3P fi=NULL); + + //e Calculates Norm matrix + //e The result is saved in AWPMD::Norm[s] + void norm_matrix(int s); + + //e Performs LU-factorization of the Norm matrix + //e AWPMD::Norm[s] is replaced by the LU matrix + void norm_factorize(int s); + + //e Invert Norm matrix + //e AWPMD::Norm[s] is replaced by the inverted matrix + void norm_invert(int s); + + //e Get the determinant of the norm-matrix for the particles with spin s + double norm_matrix_det(int s); + + //e Get the determinant logarithm of the norm-matrix for the particles with spin s + double norm_matrix_detl(int s); + + double get_energy(); + + //e makes timestep of electronic component (NOT IMPLEMENTED) + int step(double dt); + + ///\en Gets current electronic coordinates. + /// Transforms the momenta to velocity v according to mass setting (-1 means me) + int get_electrons(int spin, Vector_3P x, Vector_3P v, double* w, double* pw, double mass=-1); + + void set_harm_constr(double w0) { + constraint = HARM; + harm_w0_4 = h_sq*9./(8.*m_electron)/(w0*w0*w0*w0); + } + + void set_fix_constr(double w0_) { + constraint = FIX; + w0 = w0_; + } + + ///\en Prepares force arrays according to \a flag setting for interaction() + virtual void clear_forces(int flagi,Vector_3P fi, Vector_3P fe_x, + Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c=NULL); + + + ///\en Creates wave packet acording to the given physical parameters. + /// The function may change its arguments by applying existing constraints! + /// Default mass (-1) is the electron mass AWPMD::me. + WavePacket create_wp(Vector_3 &x, Vector_3 &v, double &w, double &pw, double mass=-1); + + +}; + + +# endif diff --git a/lib/awpmd/systems/interact/TCP/wpmd_split.cpp b/lib/awpmd/systems/interact/TCP/wpmd_split.cpp new file mode 100644 index 0000000000..361fae1f34 --- /dev/null +++ b/lib/awpmd/systems/interact/TCP/wpmd_split.cpp @@ -0,0 +1,1325 @@ +# include "wpmd_split.h" +//# include "erf.h" + + +void AWPMD_split::resize(int flag){ + for(int s=0;s<2;s++){ + wf_norm[s].resize(ne[s]); + if(flag&(0x8|0x4)){ //electron forces needed + wf_norm_der[s].resize(nvar[s]); + ovl_der[s].resize(nvar[s]); + E_der[s].resize(nvar[s]); + + if(flag&(0x8|0x4) || norm_needed){ //electron forces or norm matrix needed + if(approx==HARTREE){ // L and Norm are needed in block form + Lh[s].resize(nvar[s]); + if(norm_needed){ + Normh[s].resize(ne[s]); + for(int i=0;i 1) + return LOGERR(-1,fmt("AWPMD_split.set_electrons: invaid spin setting (%d)!",s),LINFO); + + // calculating the total n + nvar[s]=0; + int n=0; + for(int i=0;i0){ // width PBC, keeping the width are within [0,Lextra] + w[i]=fmod(w[i],Lextra); + if(w[i]<0) w[i]+=Lextra; + rw=w[i]; // WP width for energy evaluation is within [0, L/2] + if(rw > Lextra/2) rw = Lextra - rw; + } + else + rw=w[i]; + + wp[s][i].init(rw,x[i],v[i]*m_electron/h_plank,pw[i]/h_plank);*/ + wp[s][i]=create_wp(x[i],v[i],w[i],pw[i],mass); + //printf("%15d %15g\n",i,rw); + // assign default partition + partition1[s].push_back(i+1); + } + + // assign electronic charge + if(q) + qe[s].assign(q,q+nwp[s]); + else + qe[s].assign(nwp[s],-1); + + + + return 1; +} + + + +void AWPMD_split::eterm_deriv(int ic1,int s1,int c1, int j1,int ic2,int s2, int c2, int k2,cdouble pref, + const OverlapDeriv &o,cdouble v,cdouble dv_aj_conj, + cdouble dv_ak,cVector_3 dv_bj_conj, cVector_3 dv_bk){ + cdouble cj(split_c[s1][ic1+j1][0],split_c[s1][ic1+j1][1]); + cdouble ck(split_c[s2][ic2+k2][0],split_c[s2][ic2+k2][1]); + int indw1=8*ic1, indw2=8*ic2; + int indn1=(nvar[s1]/10)*8+2*ic1, indn2=(nvar[s2]/10)*8+2*ic2; + cdouble part_jk=conj(cj)*ck; + + int M= 1; //(j1==k2 ? 1 : 2); + + // over a_k_re + E_der[s2][indw2+8*k2]+= M*real(pref*part_jk*(o.da2_re()*v+o.I0*dv_ak)); + // over a_k_im + E_der[s2][indw2+8*k2+1]+= M*real(pref*part_jk*(o.da2_im()*v+i_unit*o.I0*dv_ak)); + // over a_j_re + E_der[s1][indw1+8*j1]+= M*real(pref*part_jk*(o.da1_re()*v+o.I0*dv_aj_conj)); + // over a_j_im + E_der[s1][indw1+8*j1+1]+= M*real(pref*part_jk*(o.da1_im()*v-i_unit*o.I0*dv_aj_conj)); + + for(int i=0;i<3;i++){ + // over b_k_re + E_der[s2][indw2+8*k2+2+2*i]+= M*real(pref*part_jk*(o.db2_re(i)*v+o.I0*dv_bk[i])); + // over b_k_im + E_der[s2][indw2+8*k2+2+2*i+1]+= M*real(pref*part_jk*(o.db2_im(i)*v+i_unit*o.I0*dv_bk[i])); + // over b_j_re + E_der[s1][indw1+8*j1+2+2*i]+= M*real(pref*part_jk*(o.db1_re(i)*v+o.I0*dv_bj_conj[i])); + // over b_j_im + E_der[s1][indw1+8*j1+2+2*i+1]+= M*real(pref*part_jk*(o.db1_im(i)*v-i_unit*o.I0*dv_bj_conj[i])); + } + + + // over ck_re + E_der[s2][indn2+2*k2]+=M*real(pref*conj(cj)*o.I0*v); + // over ck_im + E_der[s2][indn2+2*k2+1]+=M*real(pref*i_unit*conj(cj)*o.I0*v); + // over cj_re + E_der[s1][indn1+2*j1]+=M*real(pref*ck*o.I0*v); + // over cj_im + E_der[s1][indn1+2*j1+1]+=M*real(-pref*i_unit*ck*o.I0*v); + + double t= -M*real(pref*part_jk*o.I0*v); + // nonlocal terms: TODO: make a separate global loop for summation of nonlocal terms + for(int j=0;j -mu, v -> -v !!! + } + + + for(int k1=j1+1;k1(E_der[s1].begin()+indw1,(double *)&fe_x[iv1+k1],(double *)&fe_p[iv1+k1],&fe_w[iv1+k1],&fe_pw[iv1+k1], 1./one_h); + fe_c[iv1+k1]+=-Vector_2(E_der[s1][indn1+2*k1],E_der[s1][indn1+2*k1+1]); + }// k1 + ic1+=nspl[s1][c1]; // incrementing block1 wp address + iv1+=nspl[s1][c1]; // incrementing global variable address + } + } + } +} + +//e same as interaction, but using Hartee factorization (no antisymmetrization) +int AWPMD_split::interaction_hartree(int flag, Vector_3P fi, Vector_3P fe_x, + Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){ + + // resize arrays if needed + enum APPROX tmp=HARTREE; + swap(tmp,approx); // do not neeed large matrices + resize(flag); + swap(tmp,approx); + + // clearing forces + clear_forces(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c); + // calculate block norms and (optionally) derivatives + calc_norms(flag); + + Eee = Ew = 0.; + for(int s1=0;s1<2;s1++){ + Ee[s1]=0.; + Eei[s1]=0.; + int ic1=0; // starting index of the wp for current electron + + for(int c1=0;c1 ovl_tolerance; + } + ik+=nspl[s][k]; // incrementing block1 wp address + } + O[s] = Y[s]; // save overlap matrix + + //3. inverting the overlap matrix + int info=0; + if(nes && approx!=HARTREE){ + /*FILE *f1=fopen(fmt("matrO_%d.d",s),"wt"); + fileout(f1,Y[s],"%15g"); + fclose(f1);8*/ + + ZPPTRF("L",&nes,Y[s].arr,&info); + // analyze return code here + if(info<0) + return LOGERR(info,fmt("AWPMD.interacton: call to ZPTRF failed (exitcode %d)!",info),LINFO); + ZPPTRI("L",&nes,Y[s].arr,&info); + if(info<0) + return LOGERR(info,fmt("AWPMD.interacton: call to ZPTRI failed (exitcode %d)!",info),LINFO); + + + /*f1=fopen(fmt("matrY_%d.d",s),"wt"); + fileout(f1,Y[s],"%15g"); + fclose(f1);*/ + } + + // Clearing matrices for electronic forces + /* + if(flag&0x4){ + Te[s].Set(0.); + Tei[s].Set(0.); + }*/ + } +# if 1 + // calculating the M matrix + if(norm_needed || ( flag&(0x8|0x4) && approx!=HARTREE ) ){ + for(int s1=0;s1<2;s1++){ + if(approx!=HARTREE){ + M[s1].Set(0.); + L[s1].Set(0.); + if(norm_needed) + Norm[s1].Set(0.); + } + else + Lh[s1].assign(nvar[s1],0.); + + int indn=(nvar[s1]/10)*8; + int ic1=0; + for(int c1=0;c1, some fixes to UHF + * + * Revision 1.11 2011/05/27 08:43:52 valuev + * fixed split packet antisymmetrized version + * + * Revision 1.10 2011/05/25 05:23:43 valuev + * fixed variable transformation for norm matrix + * + * Revision 1.9 2011/05/24 19:54:32 valuev + * fixed sqmatrix::iterator + * + * Revision 1.8 2011/05/20 21:39:49 valuev + * separated norm calculation + * + * Revision 1.7 2011/05/14 18:56:19 valuev + * derivative for ee split interactions + * + * Revision 1.6 2011/05/05 08:56:02 valuev + * working split wp version + * + * Revision 1.5 2011/05/04 16:48:52 valuev + * fixed syntax + * + * Revision 1.4 2011/05/04 09:04:48 valuev + * completed wp_split (except for ee forces) + * + * Revision 1.3 2011/04/22 09:54:24 valuev + * working on split WP + * + * Revision 1.1 2011/04/20 08:43:09 valuev + * started adding split packet WPMD + * + *******************************************************************************/ + +#include "wpmd.h" + + +class AWPMD_split: public AWPMD { +protected: + int s_add, spl_add; +public: + vector split_c[2]; ///<\en split coefficients for electrons (c_re, c_im) or (psi,phi) depending on the norm mode + vector nspl[2]; ///<\en number of wave packets for each electron (size is ne[i]) + + vector wf_norm[2]; ///<\en norms for each electron + vector wf_norm_der[2]; ///<\en norm derivative + vector ovl_der[2]; ///<\en overlap derivative: \ + vector E_der[2]; ///<\en energy derivative with respect to {a,b} coordinates + + + vector< cdouble > Lh[2]; ///<\en Substitute for L in Hartree case: block matrices 1x(10*nspl[i]) + vector< sqmatrix > Normh[2]; ///<\en Substitute for Norm in Hartree case: block matrices + + ///\en resizes all internal arrays according to new electrons (wavepackets) added + virtual void resize(int flag); + + + + +public: + AWPMD_split():s_add(0),spl_add(0){} + + + ///\en Prepares to setup a new system of particles using \ref add_ion(), + /// \ref add_electron() and \ref add_split(). + /// There is no need to call this function when using + /// \ref set_electrons() and \ref set_ions() to setup particles. + virtual void reset(){ + for(int s=0;s<2;s++){ + split_c[s].clear(); + nspl[s].clear(); + } + s_add=0; + spl_add=0; + AWPMD::reset(); + } + + + + ///\en Setup electrons: forms internal wave packet representations. + /// If PBCs are used the coords must be within the range [0, cell) + /// the \a splits array defines the number of wavepackets required for each electron + /// the data for splits should be placed in the corresponding data arrays + /// \a c array contains the splits mixing coefficints + /// \a n is the number of electrons of a given spin component + /// Electron velocity v is multiplied by mass to obtain momentum. + /// Default mass (-1) means me. + /// Electronic charges q are -1 by default (when q=NULL), otherwise the charges are assigned for each split + int set_electrons(int s, int nel, Vector_3P x, Vector_3P v, double* w, double* pw, Vector_2 *c, int *splits, double mass=-1, double *q=NULL); + + + ///\en Starts adding new electron: continue with \ref add_split functions. + int add_electron(int s){ + if(s < 0 || s > 1) + return LOGERR(-1,fmt("AWPMD_split.add_electron: invaid spin setting (%d)!",s),LINFO); + s_add=s; + spl_add=0; + return ne[s_add]; + } + + ///\en Adds a new split to current electron. + /// May change the arguments according to the constraints set. + /// \return global id of the wavepacket (starting from 0 for each spin s) + /// Electron velocity v is multiplied by mass to obtain momentum. + /// Default mass (-1) means me. + /// The tags must be nonzero, >0 for the local particle, <0 for ghost particle. + /// Unique particle id is abs(tag). + /// Default tag (0) means inserting the current particle id as local particle. + int add_split(Vector_3 &x, Vector_3 &v, double &w, double &pw, Vector_2 &c, double mass=-1, double q=-1., int tag=0); + + + ///\en gets current electronic coordinates, and (optionally) number of wave packets for each electron + int get_electrons(int spin, Vector_3P x, Vector_3P v, double* w, double* pw, cdouble *c, int *splits=NULL, double mass=-1); + + + void eterm_deriv(int ic1,int s1, int c1,int k1,int ic2,int s2, int c2,int j2,cdouble pref, + const OverlapDeriv &o,cdouble v,cdouble dv_aj_conj, + cdouble dv_ak,cVector_3 dv_bj_conj, cVector_3 dv_bk); + + ///\en adds the derivatives of Y for the term v*Y[s](c2,c1) + void y_deriv(cdouble v,int s, int c2, int c1); + + + ///\en Calculates block norms an derivatives + void calc_norms(int flag); + + ///\en Prepares force arrays according to \a flag setting for interaction() + virtual void clear_forces(int flagi,Vector_3P fi, Vector_3P fe_x, + Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c); + + ///\en Calcualtes the overlap between two electrons taking all split WPs into account. + /// Norms must be pre-calculated. + cdouble overlap(int ic1, int s1, int c1,int ic2, int s2, int c2); + + //e same as interaction, but using Hartee factorization (no antisymmetrization) + int interaction_hartree(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL, + Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL); + + ///\en Calculates interaction in the system of ni ions + electrons + /// the electonic subsystem must be previously setup by set_electrons, ionic by set_ions + /// 0x1 -- give back ion forces \n + /// 0x2 -- add ion forces to the existing set \n + /// 0x4 -- calculate electronic forces \n + /// 0x8 -- add electronic forces to the existing arrays \n + /// 0x10 -- calculate internal electronic derivatives only: \n + /// will not update electronic force arrays, which may be NULL, \n + /// the forces may be obtained then using \ref get_el_forces() for all WPs \n + /// or separately for each WP using \ref get_wp_force() + /// if PBCs are used the coords must be within a range [0, cell) + int interaction(int flag=0, Vector_3P fi=NULL, Vector_3P fe_x=NULL, + Vector_3P fe_p=NULL, double *fe_w=NULL, double *fe_pw=NULL, Vector_2P fe_c=NULL); + + ///\en Get electronic forcess in the arrays provided, using calculated internal representation + /// Valid flag settings are:\n + /// 0x4 -- overwrite existing forces \n + /// 0x8 -- add electronic forces to the existing arrays \n + void get_el_forces(int flag, Vector_3P fe_x, + Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c); + + + void get_wp_force(int s, int ispl, Vector_3P fe_x, Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){ + WavePacket wk=wp[s][ispl]; + int indw1=8*ispl; + int indn1=(nvar[s]/10)*8+2*ispl; + wk.int2phys_der< eq_minus_second >(E_der[s].begin()+indw1,(double *)fe_x,(double *)fe_p,fe_w,fe_pw,1./one_h); + *fe_c=-Vector_2(E_der[s][indn1],E_der[s][indn1+1]); + } +}; + + + + + +# endif