diff --git a/CMakeLists.txt b/CMakeLists.txt index cc37429ff..408d4ec6e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -167,8 +167,7 @@ endif() set (MMAP_LIB lib/src/readlog.cpp lib/src/testtraiettoria.cpp - lib/src/traiettoria.cpp - lib/src/greenkubo2componentionicfluid.cpp + lib/src/trajectory.cpp lib/src/posizioniequilibrio.cpp lib/src/chargefluxts.cpp ) @@ -177,7 +176,6 @@ set (ANALISI_LIB lib/src/calcoliblocchi.cpp lib/src/convertibinario.cpp lib/src/convolution.cpp - lib/src/correlatorespaziale.cpp lib/src/cronometro.C lib/src/doubleround.cpp lib/src/gofrt.cpp @@ -190,15 +188,12 @@ set (ANALISI_LIB lib/src/msd.cpp lib/src/rnd.c lib/src/spettrovibrazionale.cpp - lib/src/heatc.cpp lib/src/calcolamultithread.cpp - lib/src/centerdiff.cpp - lib/src/centerofmassdiff.cpp lib/src/specialfunctions.cpp lib/src/sphericalcorrelations.cpp lib/src/steinhardt.cpp lib/src/atomicdensity.cpp - lib/src/traiettoriabase.cpp + lib/src/basetrajectory.cpp lib/src/neighbour.cpp lib/src/sphericalbase.cpp ) @@ -219,7 +214,7 @@ find_package(PythonLibsNew ${PYBIND11_PYTHON_VERSION} REQUIRED) if (${PYTHONLIBS_FOUND}) set (ANALISI_LIB "${ANALISI_LIB}" - lib/src/traiettoria_numpy.cpp + lib/src/trajectory_numpy.cpp lib/src/readlog_numpy.cpp ) endif() endif() diff --git a/analisi/main.cpp b/analisi/main.cpp index 38be90d67..a4c01a28b 100644 --- a/analisi/main.cpp +++ b/analisi/main.cpp @@ -22,7 +22,6 @@ #include #include "config.h" #include "istogrammavelocita.h" -#include "greenkubo2componentionicfluid.h" #include "greenkuboNcomponentionicfluid.h" #include "convolution.h" #include "heatfluxts.h" @@ -49,6 +48,8 @@ #endif #include +#include "blockaverage.h" + namespace std{ template @@ -160,7 +161,7 @@ int main(int argc, char ** argv) int sub_mean_start=0,numero_frame=0,every=1,blocksize=0,elast=0,blocknumber=0,numero_thread,nbins_vel,skip=1,conv_n=20,final=60,stop_acf=0; unsigned int n_seg=0,gofrt=0,read_lines_thread=200,sph=0,buffer_size=10; bool sub_mean=false,test=false,spettro_vibraz=false,velocity_h=false,heat_coeff=false,debug=false,debug2=false,dumpGK=false,msd=false,msd_cm=false,msd_self=false,bench=false,fpe=false; - double vmax_h=0,cariche[2],dt=5e-3,vicini_r=0.0; + double vmax_h=0,charge[2],dt=5e-3,vicini_r=0.0; std::pair nk; std::array kdir; std::vector cvar_list,kk_l; @@ -187,8 +188,8 @@ int main(int argc, char ** argv) ("skip,s",boost::program_options::value(&skip)->default_value(1),"when an average over the trajectory is performed, consecutive steps have a distance specified with this option") ("every,e",boost::program_options::value(&every)->default_value(1),"For sums over time lags every steps have a difference multiple of every. So far implemented only for g(r,t) and give the number of steps between two consecutive t.") #ifdef EXPERIMENTAL - ("charge1",boost::program_options::value(&cariche[0])->default_value(1.0),"charge of type 1 (used in a particular case)") - ("charge2",boost::program_options::value(&cariche[1])->default_value(-1.0),"charge of type 2 (used in a particular case)") + ("charge1",boost::program_options::value(&charge[0])->default_value(1.0),"charge of type 1 (used in a particular case)") + ("charge2",boost::program_options::value(&charge[1])->default_value(-1.0),"charge of type 2 (used in a particular case)") ("conv_n,C",boost::program_options::value(&conv_n)->default_value(10),"sigma of the gaussian that may be used to compute a convolution with the green-kubo integrals, to smooth out some noise. Number of points units.") ("final,f",boost::program_options::value(&final)->default_value(60),"number of points to use to extract the final value of gk integral (used in a particular case)") #endif @@ -301,9 +302,9 @@ int main(int argc, char ** argv) if (debug2){ ReadLog<> test(log_input,0,1,numero_thread,read_lines_thread,headers); - Traiettoria * binary_traj=NULL; + Trajectory * binary_traj=NULL; if (test.need_binary(headers)>0) { - binary_traj=new Traiettoria(input); + binary_traj=new Trajectory(input); test.calc_currents(binary_traj,blocknumber); } for (unsigned int i=0;i test(log_input,0,1,numero_thread,read_lines_thread,headers); - Traiettoria * binary_traj=NULL; + Trajectory * binary_traj=NULL; //qui devo aggiungere la traiettoria binaria a ReadLog, qualora ReadLog ne constati la necessità if (test.need_binary(headers)>0) { - binary_traj=new Traiettoria(input); + binary_traj=new Trajectory(input); test.calc_currents(binary_traj,blocknumber); } @@ -390,68 +391,12 @@ int main(int argc, char ** argv) if (headers.size()==0){ - - MediaBlocchiG,GreenKubo2ComponentIonicFluid,std::string,double*,unsigned int,bool,unsigned int,unsigned int> - greenK_c(&test,blocknumber); - - MediaVarCovar greenK(GreenKubo2ComponentIonicFluid::narr,cvar); - - greenK_c.calcola_custom(&greenK,log_input,cariche,skip,dumpGK,stop_acf,numero_thread); - - double factors[GreenKubo2ComponentIonicFluid::narr]={ - factor_conv*factor_intToCorr, //Jee - factor_conv2*factor_intToCorr, //Jzz - factor_conv*factor_intToCorr, //Jez - factor_conv, //intee - factor_conv2, //intzz - factor_conv, //intez - factor_conv, //lambda - factor_conv*factor_intToCorr, //Jze - factor_conv, //intze - factor_conv, //int_ein_ee - factor_conv, //int_ein_ez - factor_conv, //int_ein_ze - factor_conv2, //int_ein_zz - factor_conv //lambda_einst - }; - double *lambda_conv=new double[greenK.size()]; - double *lambda_conv_var=new double[greenK.size()]; - - if (final>=greenK.size()) final=greenK.size()-1; - if (final <0) final=0; - - Convolution convoluzione( std::function ([&conv_n](const double & x)->double{ - return exp(-x*x/(2*conv_n*conv_n)); - }),(conv_n*6+1),-3*conv_n,3*conv_n,3*conv_n); - convoluzione.calcola(greenK.media(6),lambda_conv,greenK.size(),1); - convoluzione.calcola(greenK.varianza(6),lambda_conv_var,greenK.size(),1); - if (factors_input.size()==0) - std::cout << "# T T1sigma factor1 factor2\n#" - <,GreenKuboNComponentIonicFluid >, + BlockAverageG,GreenKuboNComponentIonicFluid >, std::string, unsigned int, std::vector, @@ -530,11 +475,11 @@ int main(int argc, char ** argv) } } - Traiettoria test(input); + Trajectory test(input); #define MSD_(fpe_)\ - {using MSD=MSD;\ - MediaBlocchi Msd(&test,blocknumber);\ - Msd.calcola(skip,stop_acf,numero_thread,msd_cm,msd_self,dumpGK);\ + {using MSD=MSD;\ + BlockAverage Msd(&test,blocknumber);\ + Msd.calculate(skip,stop_acf,numero_thread,msd_cm,msd_self,dumpGK);\ for (unsigned int i=0;ilunghezza()/test.get_ntypes()/f_cm;i++) {\ for (unsigned int j=0;jelemento(i*test.get_ntypes()*f_cm+j) << " " <<\ @@ -552,12 +497,12 @@ int main(int argc, char ** argv) throw std::runtime_error("You have to specify the distance range with the option -F.\n"); } std::cerr << "Calculation of g(r,t) -- distinctive and non distinctive part of the van Hove function...\n"; - Traiettoria tr(input); + Trajectory tr(input); tr.set_pbc_wrap(true); //è necessario impostare le pbc per far funzionare correttamente la distanza delle minime immagini - MediaBlocchi,double,double,unsigned int,unsigned int,unsigned int, unsigned int,unsigned int, bool> + BlockAverage,double,double,unsigned int,unsigned int,unsigned int, unsigned int,unsigned int, bool> gofr(&tr,blocknumber); - gofr.calcola(factors_input[0],factors_input[1],gofrt,stop_acf,numero_thread,skip,every,dumpGK); + gofr.calculate(factors_input[0],factors_input[1],gofrt,stop_acf,numero_thread,skip,every,dumpGK); unsigned int ntyp=tr.get_ntypes()*(tr.get_ntypes()+1); unsigned int tmax=gofr.media()->lunghezza()/gofrt/ntyp; @@ -587,20 +532,20 @@ int main(int argc, char ** argv) throw std::runtime_error("You have to specify the distance range with the option -F.\n"); } std::cerr << "Calculation of spherical harmonic density correlation function is beginning (this can take a lot of time)...\n"; - Traiettoria tr(input); + Trajectory tr(input); tr.set_pbc_wrap(false); //è necessario impostare le pbc per far funzionare correttamente la distanza delle minime immagini - using SHC=SphericalCorrelations<10,double,Traiettoria>; + using SHC=SphericalCorrelations<10,double,Trajectory>; SHC::rminmax_t rminmax; for (unsigned int i=0; i sh(&tr,blocknumber); - sh.calcola(rminmax,sph,stop_acf,numero_thread,skip,buffer_size,dumpGK,{}); + sh.calculate(rminmax,sph,stop_acf,numero_thread,skip,buffer_size,dumpGK,{}); auto shape= sh.media()->get_shape(); @@ -617,18 +562,18 @@ int main(int argc, char ** argv) }else if (vicini_r>0){ std::cerr << "Beginning of calculation of neighbour histogram\n"; - Traiettoria test(input); + Trajectory test(input); test.set_pbc_wrap(true); IstogrammaAtomiRaggio h(&test,vicini_r,skip,numero_thread); unsigned int nt=test.get_ntimesteps(), s=(nt-1)/blocknumber; h.reset(s); - test.imposta_dimensione_finestra_accesso(s); + test.set_data_access_block_size(s); for (unsigned int i=0;i *hist=h.get_hist(); for (unsigned int i=0;i0) { std::cerr << "Velocity histogram...\n"; - Traiettoria test(input); - MediaBlocchi istogramma_vel(&test,blocknumber); - istogramma_vel.calcola(nbins_vel,vmax_h); + Trajectory test(input); + BlockAverage istogramma_vel(&test,blocknumber); + istogramma_vel.calculate(nbins_vel,vmax_h); //stampa i risultati unsigned int lungh_sing_h=istogramma_vel.media()->lunghezza()/(3*test.get_ntypes()); @@ -658,20 +603,20 @@ int main(int argc, char ** argv) } else if (fononefile != "") { std::cerr << "Normal mode coordinate transformation...\n"; - Traiettoria test(input); + Trajectory test(input); blocksize = test.get_ntimesteps()/blocknumber; ModiVibrazionali test_normali(&test,ifcfile,fononefile,numero_thread,blocksize); test_normali.reset(stop_acf); - test_normali.calcola(0); + test_normali.calculate(0); return 1; } else if (spettro_vibraz) { std::cerr << "Vibrational spectrum...\n"; - Traiettoria test(input); + Trajectory test(input); // SpettroVibrazionale test_spettro(&test); - MediaBlocchi,bool> test_spettro_blocchi(&test,blocknumber); + BlockAverage,bool> test_spettro_blocchi(&test,blocknumber); - test_spettro_blocchi.calcola(dumpGK); + test_spettro_blocchi.calculate(dumpGK); for (unsigned int i=0;ilunghezza()/(3*test.get_ntypes());i++) { std::cout << i << " "; for (unsigned int itipo=0;itipo #include template -class AtomicDensity : public CalcolaMultiThread>, public OperazioniSuLista,Hist> +class AtomicDensity : public CalculateMultiThread>, public VectorOp,Hist> { public: AtomicDensity(T *t, @@ -18,12 +18,12 @@ class AtomicDensity : public CalcolaMultiThread>, public O ~AtomicDensity(); void calc_single_th(const unsigned int &start, const unsigned int & stop, const unsigned int & primo, const unsigned int & ith); void join_data(); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + unsigned int nExtraTimesteps(unsigned int n_b); void reset(const unsigned int numeroTimestepsPerBlocco) ; std::vector get_shape() const; std::vector get_stride() const; using This = AtomicDensity; - using OperazioniSuLista::azzera; + using VectorOp::azzera; private: @@ -52,11 +52,11 @@ class AtomicDensity : public CalcolaMultiThread>, public O return idx(idxs[2],idxs[1],idxs[0],itype); } - using CalcolaMultiThread::ntimesteps; - using CalcolaMultiThread::skip; - using CalcolaMultiThread::nthreads; - using OperazioniSuLista::lunghezza_lista; - using OperazioniSuLista::lista; + using CalculateMultiThread::ntimesteps; + using CalculateMultiThread::skip; + using CalculateMultiThread::nthreads; + using VectorOp::data_length; + using VectorOp::vdata; }; diff --git a/lib/include/traiettoriabase.h b/lib/include/basetrajectory.h similarity index 67% rename from lib/include/traiettoriabase.h rename to lib/include/basetrajectory.h index 4cc7e9c7c..5584d8e7c 100644 --- a/lib/include/traiettoriabase.h +++ b/lib/include/basetrajectory.h @@ -43,26 +43,26 @@ struct Cell : public Cell_ { }; template -class TraiettoriaBase { +class BaseTrajectory { public: enum class BoxFormat{Lammps_ortho,Lammps_triclinic,Cell_vectors, Invalid}; - TraiettoriaBase() : ntypes{0},n_timesteps{0},natoms{0},min_type{0},max_type{0},timestep_corrente{0}, - wrap_pbc{true}, buffer_posizioni{nullptr}, buffer_velocita{nullptr}, - buffer_scatola{nullptr}, buffer_posizioni_cm{nullptr}, - buffer_velocita_cm{nullptr}, masse{nullptr}, cariche{nullptr}, - buffer_tipi{nullptr},buffer_tipi_id{nullptr}, serve_pos{true}, + BaseTrajectory() : ntypes{0},n_timesteps{0},natoms{0},min_type{0},max_type{0},current_timestep{0}, + wrap_pbc{true}, buffer_positions{nullptr}, buffer_velocity{nullptr}, + buffer_boxes{nullptr}, buffer_positions_cm{nullptr}, + buffer_velocity_cm{nullptr}, mass{nullptr}, charge{nullptr}, + buffer_type{nullptr},buffer_type_id{nullptr}, serve_pos{true}, box_format{BoxFormat::Invalid}{} - //double * posizioni (const int & timestep, const int & atomo) { return static_cast(this)->posizioni(timestep,atomo);} - DECL_CALL_BASE_2(double *, posizioni, (const int &, timestep), (const int &, atomo)) - DECL_CALL_BASE_2(double *, velocita,(const int &, timestep), (const int &, atomo) ) - DECL_CALL_BASE_1(double *, scatola, (const int &, timestep)) - DECL_CALL_BASE_2(double *, posizioni_cm,(const int &, timestep), (const int &, tipo)) - DECL_CALL_BASE_2(double *, velocita_cm,(const int &, timestep), (const int &, tipo)) - DECL_CALL_BASE_0(double *,scatola_last) + //double * positions (const int & timestep, const int & atomo) { return static_cast(this)->positions(timestep,atomo);} + DECL_CALL_BASE_2(double *, positions, (const int &, timestep), (const int &, atomo)) + DECL_CALL_BASE_2(double *, velocity,(const int &, timestep), (const int &, atomo) ) + DECL_CALL_BASE_1(double *, box, (const int &, timestep)) + DECL_CALL_BASE_2(double *, positions_cm,(const int &, timestep), (const int &, tipo)) + DECL_CALL_BASE_2(double *, velocity_cm,(const int &, timestep), (const int &, tipo)) + DECL_CALL_BASE_0(double *,box_last) std::vector get_types(){ get_ntypes(); return types; @@ -70,15 +70,14 @@ class TraiettoriaBase { } unsigned int get_type(const unsigned int &atomo){ if (atomo < natoms) { - return buffer_tipi_id[atomo]; + return buffer_type_id[atomo]; } else { throw std::runtime_error("Atom index out of range\n"); } } enum Errori {non_inizializzato=0,oltre_fine_file=2,Ok=1}; - Errori imposta_dimensione_finestra_accesso(const size_t & timesteps){std::cerr << "Warning: doing nothing (not reading in blocks)"< void pbc_wrap(ssize_t idx) { - double * c = buffer_scatola+buffer_scatola_stride*idx; + double * c = buffer_boxes+buffer_boxes_stride*idx; Cell s(c); double x0[3]; s.middle(x0); for (size_t iatom=0;iatom(itimestep,i); - double *xj=posizioni(jtimestep,j); - const double *l=scatola(itimestep)+3; - const double *xy_xz_yz = scatola(itimestep)+6; + double *xi=positions(itimestep,i); + double *xj=positions(jtimestep,j); + const double *l=box(itimestep)+3; + const double *xy_xz_yz = box(itimestep)+6; d2=d2_minImage_triclinic(xi,xj,l,x,xy_xz_yz); return d2; } + /** + * This returns the square of the minimum image distance between the specified timesteps and atom + * using the box data specified in the arguments (if not triclinc the last argument is not used) + * Sets x(i) - x(j) in memory located at x + * **/ double d2_minImage_triclinic(double * xi, ///position of first atom double *xj, ///position of second atom const double *l, /// half orthorombic edges @@ -194,6 +217,10 @@ class TraiettoriaBase { } return d2; } + /** + * apply the minimum image convention in the triclinic and not triclinc case. + * Must provide box data as in the internal format used by the code. + * **/ template static void minImage_triclinic(double * __restrict delta, const double * __restrict l_half, @@ -241,7 +268,7 @@ class TraiettoriaBase { } - size_t get_size_posizioni() const { return buffer_posizioni_size; } + size_t get_positions_size() const { return buffer_positions_size; } size_t get_size_cm() const {return buffer_cm_size;} std::vector get_shape() { return {static_cast(loaded_timesteps), @@ -270,7 +297,7 @@ class TraiettoriaBase { return triclinic; } - size_t get_box_stride() const {return buffer_scatola_stride;} + size_t get_box_stride() const {return buffer_boxes_stride;} //functions that are compiled in cpp file void dump_lammps_bin_traj(const std::string &fname, int start_ts, int stop_ts); @@ -280,19 +307,19 @@ class TraiettoriaBase { protected: - double * buffer_posizioni; //velocita' e posizioni copiate dal file caricato con mmap, in ordine (nela traiettoria di LAMMPS sono disordinate) - double * buffer_velocita; - double * masse; - double * cariche; - double * buffer_scatola; //dimensioni della simulazione ad ogni timestep - size_t buffer_scatola_stride; // 6 for orthorombic, 9 for triclinic + double * buffer_positions; //positions and velocities may need a copy because in the lammps format they are not ordered + double * buffer_velocity; + double * mass; //used almost nowhere + double * charge; + double * buffer_boxes; //box data in specified format + size_t buffer_boxes_stride; // 6 for orthorombic, 9 for triclinic BoxFormat box_format; //format used to store box information - double * buffer_posizioni_cm; // posizioni del centro di massa - double * buffer_velocita_cm; // velocità del centro di massa - size_t buffer_posizioni_size, buffer_cm_size; //sizes of allocated buffers + double * buffer_positions_cm; // center of mass position for each atomic type + double * buffer_velocity_cm; // center of mass velocity for each atomic type + size_t buffer_positions_size, buffer_cm_size; //sizes of allocated buffers - int * buffer_tipi,*buffer_tipi_id; - ssize_t natoms,ntypes,min_type,max_type,n_timesteps, loaded_timesteps,timestep_corrente; + int * buffer_type,*buffer_type_id; + ssize_t natoms,ntypes,min_type,max_type,n_timesteps, loaded_timesteps,current_timestep; bool wrap_pbc, triclinic; diff --git a/lib/include/mediablocchi.h b/lib/include/blockaverage.h similarity index 71% rename from lib/include/mediablocchi.h rename to lib/include/blockaverage.h index b0bd07b34..19fe4fcf7 100644 --- a/lib/include/mediablocchi.h +++ b/lib/include/blockaverage.h @@ -13,7 +13,7 @@ #ifndef MEDIABLOCCHI_H #define MEDIABLOCCHI_H -#include "traiettoria.h" +#include "trajectory.h" #include "assert.h" #include #include @@ -31,20 +31,20 @@ template class TraiettoriaF { public: - static void imposta_dimensione_finestra_accesso(unsigned int s_,TR* traiettoria ){} - static void imposta_inizio_accesso(unsigned int s_,TR* traiettoria ){} + static void set_data_access_block_size(unsigned int s_,TR* traiettoria ){} + static void set_access_at(unsigned int s_,TR* traiettoria ){} static unsigned int get_ntimesteps(TR* traiettoria ){abort();return 0;} }; -template <> class TraiettoriaF { +template <> class TraiettoriaF { public: -static void imposta_dimensione_finestra_accesso(unsigned int s_,Traiettoria* traiettoria ){ - traiettoria->imposta_dimensione_finestra_accesso(s_); +static void set_data_access_block_size(unsigned int s_,Trajectory* traiettoria ){ + traiettoria->set_data_access_block_size(s_); } -static void imposta_inizio_accesso(unsigned int s_,Traiettoria* traiettoria ){ - traiettoria->imposta_inizio_accesso(s_); +static void set_access_at(unsigned int s_,Trajectory* traiettoria ){ + traiettoria->set_access_at(s_); } -static unsigned int get_ntimesteps(Traiettoria* traiettoria ){ +static unsigned int get_ntimesteps(Trajectory* traiettoria ){ return traiettoria->get_ntimesteps(); } @@ -53,10 +53,10 @@ static unsigned int get_ntimesteps(Traiettoria* traiettoria ){ template <> class TraiettoriaF > { public: -static void imposta_dimensione_finestra_accesso(unsigned int s_,ReadLog<>* traiettoria ){ +static void set_data_access_block_size(unsigned int s_,ReadLog<>* traiettoria ){ } -static void imposta_inizio_accesso(unsigned int s_,ReadLog<>* traiettoria ){ +static void set_access_at(unsigned int s_,ReadLog<>* traiettoria ){ } static unsigned int get_ntimesteps(ReadLog<>* traiettoria ){ @@ -67,10 +67,10 @@ static unsigned int get_ntimesteps(ReadLog<>* traiettoria ){ template <> class TraiettoriaF > { public: -static void imposta_dimensione_finestra_accesso(unsigned int s_,ReadLog* traiettoria ){ +static void set_data_access_block_size(unsigned int s_,ReadLog* traiettoria ){ } -static void imposta_inizio_accesso(unsigned int s_,ReadLog* traiettoria ){ +static void set_access_at(unsigned int s_,ReadLog* traiettoria ){ } static unsigned int get_ntimesteps(ReadLog* traiettoria ){ @@ -80,13 +80,12 @@ static unsigned int get_ntimesteps(ReadLog* traiettoria ){ }; -template class MediaBlocchiG +template class BlockAverageG { public: - MediaBlocchiG(TR * t, + BlockAverageG(TR * t, const unsigned int & numero_blocchi ) { -// static_assert(std::is_base_of::value,"T deve essere derivato da Calcolo!"); traiettoria=t; n_b=numero_blocchi; ok=false; @@ -98,10 +97,7 @@ template class MediaBlocchiG } - ~MediaBlocchiG() { -#ifdef DEBUG - std::cerr << "~MediaBlocchiG(): Tmedio="< class MediaBlocchiG calcolo = new T (traiettoria,arg...); - int timestepsPerBlocco=(TraiettoriaF::get_ntimesteps(traiettoria)-calcolo->numeroTimestepsOltreFineBlocco(n_b))/n_b; + int timestepsPerBlocco=(TraiettoriaF::get_ntimesteps(traiettoria)-calcolo->nExtraTimesteps(n_b))/n_b; if(timestepsPerBlocco>0){ s=timestepsPerBlocco; ok=true; @@ -124,7 +120,7 @@ template class MediaBlocchiG calcolo->reset(s); calc->calcola_begin(s,calcolo); - TraiettoriaF::imposta_dimensione_finestra_accesso(s+calcolo->numeroTimestepsOltreFineBlocco(n_b),traiettoria); + TraiettoriaF::set_data_access_block_size(s+calcolo->nExtraTimesteps(n_b),traiettoria); cronometro cron; #ifndef USE_MPI cron.set_expected(1.0/double(n_b)); @@ -136,10 +132,10 @@ template class MediaBlocchiG #endif calcolo->reset(s); - TraiettoriaF::imposta_inizio_accesso(iblock*s,traiettoria); - calcolo->calcola(iblock*s); + TraiettoriaF::set_access_at(iblock*s,traiettoria); + calcolo->calculate(iblock*s); - calc->calcola(calcolo); + calc->calculate(calcolo); cron.stop(); #ifdef DEBUG std::cerr << "Time for block "< class MediaBlocchiG break; } calcolo->reset(s); - TraiettoriaF::imposta_inizio_accesso(ib*s,traiettoria); - calcolo->calcola(ib*s); + TraiettoriaF::set_access_at(ib*s,traiettoria); + calcolo->calculate(ib*s); cronometro cronmpi; unsigned int mpirecv=1; if (mpime==0){ //ricevi i risultati degli altri e calcola ciò che è da calcolare. - calc->calcola(calcolo); + calc->calculate(calcolo); if (iblock+mpisize class MediaBlocchiG cronmpi.start(); Mp::mpi().recv_root(calcolo,i); cronmpi.stop(); - calc->calcola(calcolo); + calc->calculate(calcolo); } } else { Mp::mpi().send_to_root(calcolo); @@ -195,7 +191,7 @@ template class MediaBlocchiG - void calcola(Args ... arg) { + void calculate(Args ... arg) { Tmedio = new T(traiettoria, arg...); Tvar = new T(traiettoria,arg...); delta = new T(traiettoria,arg...); @@ -222,7 +218,7 @@ template class MediaBlocchiG }; -template using MediaBlocchi = MediaBlocchiG; +template using BlockAverage = BlockAverageG; diff --git a/lib/include/calcoliblocchi.h b/lib/include/calcoliblocchi.h index 188d1310b..20bd4aa4b 100644 --- a/lib/include/calcoliblocchi.h +++ b/lib/include/calcoliblocchi.h @@ -32,7 +32,7 @@ template class MediaVar{ iblock=0; } - void calcola(T* calcolo) { + void calculate(T* calcolo) { #ifdef DEBUG2 std::cerr << "*delta = *calcolo - *Tmedio;\n"; #endif @@ -133,7 +133,7 @@ template class MediaVarCovar{ delete [] dx2; } - void calcola(T* calcolo) { + void calculate(T* calcolo) { assert(dx2!=0 && dx!=0); assert(calcolo->lunghezza()%n_columns==0); assert(size_block==calcolo->lunghezza()/n_columns); diff --git a/lib/include/calcolamultithread.h b/lib/include/calculatemultithread.h similarity index 77% rename from lib/include/calcolamultithread.h rename to lib/include/calculatemultithread.h index 5943e3ce2..a9fd712a8 100644 --- a/lib/include/calcolamultithread.h +++ b/lib/include/calculatemultithread.h @@ -15,7 +15,7 @@ #include -namespace CalcolaMultiThread_Flags { +namespace CalculateMultiThread_Flags { //mutually exclusive SPLITS: constexpr int PARALLEL_SPLIT_AVERAGE = 0b00000001; // assign to the single thread worker different averages range constexpr int PARALLEL_SPLIT_TIME = 0b00000010; // assigne to the single thread worker different time differences ranges @@ -34,15 +34,15 @@ constexpr int CALL_CALC_INIT = 0b100000000; } -template -class CalcolaMultiThread +template +class CalculateMultiThread { public: - CalcolaMultiThread(const ssize_t nthreads=0, const ssize_t skip=0, const size_t natoms=0, const ssize_t every=0) : nthreads{nthreads},skip{skip},ntimesteps{0},natoms{natoms},every{every} + CalculateMultiThread(const ssize_t nthreads=0, const ssize_t skip=0, const size_t natoms=0, const ssize_t every=0) : nthreads{nthreads},skip{skip},ntimesteps{0},natoms{natoms},every{every} { - if (nthreads==0) CalcolaMultiThread::nthreads=1; - if (skip==0) CalcolaMultiThread::skip=1; - if (every==0) CalcolaMultiThread::every=1 ; + if (nthreads==0) CalculateMultiThread::nthreads=1; + if (skip==0) CalculateMultiThread::skip=1; + if (every==0) CalculateMultiThread::every=1 ; } static constexpr int FLAGS=FLAGS_T; @@ -52,24 +52,24 @@ class CalcolaMultiThread res.first=ith*npassith; if (ith==nthreads-1) { res.second=end; - if constexpr ( !!(FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_AVERAGE)){ + if constexpr ( !!(FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_AVERAGE)){ //align the last timestep so it is a multiple of skip and the total number of points is ntimestep/skip res.second = end - end%skip; } - if constexpr ( !!(FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_TIME)) { + if constexpr ( !!(FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_TIME)) { //align the last timestep so it is a multiple of every and the total number of points is leff/every res.second= end - end%every; } } else { res.second=(ith+1)*npassith; } - if constexpr ( !!(FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_AVERAGE)) { + if constexpr ( !!(FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_AVERAGE)) { if (skip > 1 && res.first % skip > 0) { //align the first timestep so it is a multiple of skip res.first = res.first - res.first % skip + skip; } } - if constexpr ( !!(FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_TIME)) { + if constexpr ( !!(FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_TIME)) { if (every > 1 && res.first % every > 0 ) { res.first = res.first - res.first % every + every ; } @@ -78,34 +78,34 @@ class CalcolaMultiThread } void init_split() { - if constexpr ( !! (FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_AVERAGE)){ + if constexpr ( !! (FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_AVERAGE)){ npassith=ntimesteps/nthreads; end=ntimesteps; - } else if (FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_TIME) { + } else if (FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_TIME) { npassith=leff/nthreads; end=leff; - } else if (FLAGS & CalcolaMultiThread_Flags::PARALLEL_SPLIT_ATOM) { + } else if (FLAGS & CalculateMultiThread_Flags::PARALLEL_SPLIT_ATOM) { npassith=natoms/nthreads; end=natoms; } else { - static_assert (FLAGS & (CalcolaMultiThread_Flags::PARALLEL_SPLIT_ATOM | CalcolaMultiThread_Flags::PARALLEL_SPLIT_TIME | CalcolaMultiThread_Flags::PARALLEL_SPLIT_AVERAGE), - "YOU MUST SPECIFY ONE OF CalcolaMultiThread_Flags" ); + static_assert (FLAGS & (CalculateMultiThread_Flags::PARALLEL_SPLIT_ATOM | CalculateMultiThread_Flags::PARALLEL_SPLIT_TIME | CalculateMultiThread_Flags::PARALLEL_SPLIT_AVERAGE), + "YOU MUST SPECIFY ONE OF CalculateMultiThread_Flags" ); } t0=0;t1=1; i0=0;i1=1; - if constexpr ( !!(FLAGS & CalcolaMultiThread_Flags::SERIAL_LOOP_AVERAGE)) { + if constexpr ( !!(FLAGS & CalculateMultiThread_Flags::SERIAL_LOOP_AVERAGE)) { i0=0; i1=ntimesteps; } - if constexpr ( !!(FLAGS & CalcolaMultiThread_Flags::SERIAL_LOOP_TIME)) { + if constexpr ( !!(FLAGS & CalculateMultiThread_Flags::SERIAL_LOOP_TIME)) { t0=0; t1=leff; } } - void calcola(size_t primo){ + void calculate(size_t primo){ init_split(); - if constexpr (!!(FLAGS & CalcolaMultiThread_Flags::CALL_CALC_INIT)) { + if constexpr (!!(FLAGS & CalculateMultiThread_Flags::CALL_CALC_INIT)) { static_cast(this)->calc_init(primo); } std::exception_ptr * thread_exception = new std::exception_ptr[nthreads]; @@ -120,11 +120,11 @@ class CalcolaMultiThread // calculate given start and stop timestep. note that & captures everything, user must take care of multithread safety of calc_single_th function // select a different signature of the calculation function depending on where the loops (if present) are parallelized // rangeA and rangeB are passed to each thread and their meaning depends on what PARALLEL_SPLIT_* option you choose - if constexpr ( !! (FLAGS & CalcolaMultiThread_Flags::SERIAL_LOOP_AVERAGE && (FLAGS & CalcolaMultiThread_Flags::SERIAL_LOOP_TIME))) { + if constexpr ( !! (FLAGS & CalculateMultiThread_Flags::SERIAL_LOOP_AVERAGE && (FLAGS & CalculateMultiThread_Flags::SERIAL_LOOP_TIME))) { static_cast(this)->calc_single_th(t,i,range.first,range.second,primo,ith); // time, average, rangeA, rangeB, primo, thread id - } else if constexpr( !! (FLAGS & CalcolaMultiThread_Flags::SERIAL_LOOP_AVERAGE) ) { + } else if constexpr( !! (FLAGS & CalculateMultiThread_Flags::SERIAL_LOOP_AVERAGE) ) { static_cast(this)->calc_single_th(i,range.first,range.second,primo,ith); // average, rangeA, rangeB, primo, thread id - } else if constexpr( !! (FLAGS & CalcolaMultiThread_Flags::SERIAL_LOOP_TIME) ) { + } else if constexpr( !! (FLAGS & CalculateMultiThread_Flags::SERIAL_LOOP_TIME) ) { static_cast(this)->calc_single_th(t,range.first,range.second,primo,ith); // time, rangeA, rangeB, primo, thread id } else { static_cast(this)->calc_single_th(range.first,range.second,primo,ith); // rangeA, rangeB, primo, thread id @@ -149,12 +149,12 @@ class CalcolaMultiThread } } threads.clear(); - if constexpr (!!(FLAGS & CalcolaMultiThread_Flags::CALL_INNER_JOIN_DATA)) + if constexpr (!!(FLAGS & CalculateMultiThread_Flags::CALL_INNER_JOIN_DATA)) static_cast(this)->join_data(); } } - if constexpr (!!(FLAGS & CalcolaMultiThread_Flags::CALL_DEBUG_ROUTINE)) { + if constexpr (!!(FLAGS & CalculateMultiThread_Flags::CALL_DEBUG_ROUTINE)) { static_cast(this)->calc_end(); } delete [] thread_exception; @@ -181,7 +181,7 @@ class CalcolaMultiThread /* * those must be implemented to interface with block averages stuff. Anyway, there are useful */ - //unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b)=0; + //unsigned int nExtraTimesteps(unsigned int n_b)=0; /* * This function must allocate all the memory: The argument is the number of timesteps to be calculated. * remember to set ntimesteps here, so the division of the work can be done correctly diff --git a/lib/include/centerdiff.h b/lib/include/centerdiff.h deleted file mode 100644 index 59f9e3f12..000000000 --- a/lib/include/centerdiff.h +++ /dev/null @@ -1,30 +0,0 @@ -#ifndef CENTERDIFF_H -#define CENTERDIFF_H - -#include "calcolamultithread.h" -#include "operazionisulista.h" -#include - -class CenterDiff : public CalcolaMultiThread, public OperazioniSuLista -{ -public: - CenterDiff(Traiettoria *t, unsigned int nthreads=0, unsigned int skip=1, unsigned int nit=1,bool sum_first_two_and_ignore_vz=false,bool sum_1=false); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b) {return 0;} - void reset(const unsigned int numeroTimestepsPerBlocco); - void calc_single_th(const unsigned int &start, const unsigned int &stop, const unsigned int &primo, const unsigned int & ith) noexcept; - std::vector get_shape() const {return {lunghezza_lista/3/nit/3,nit,3,3}; } - std::vector get_stride() const {return { static_cast (sizeof (double)*3*3*nit),sizeof (double)*3*3,sizeof (double)*3,sizeof (double)};} - void set_starting_center(const std::array & s) {starting_center=s;} - void join_data(){} - ~CenterDiff(); -private: - Traiettoria *t; - using OperazioniSuLista::lista; - using OperazioniSuLista::lunghezza_lista; - unsigned int nit,lista_alloc; - bool sum_first_two_and_ignore_vz,sum_1; - std::array starting_center; - -}; - -#endif // CENTERDIFF_H diff --git a/lib/include/centerofmassdiff.h b/lib/include/centerofmassdiff.h deleted file mode 100644 index a1dea4969..000000000 --- a/lib/include/centerofmassdiff.h +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef CENTEROFMASSDIFF_H -#define CENTEROFMASSDIFF_H - -#include "calcolamultithread.h" -#include "operazionisulista.h" -#include - -class CenterOfMassDiff : public CalcolaMultiThread, public OperazioniSuLista -{ -public: - CenterOfMassDiff(Traiettoria *t, unsigned int nthreads=0, unsigned int skip=1); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b) {return 0;} - void reset(const unsigned int numeroTimestepsPerBlocco); - void calc_single_th(const unsigned int &start, const unsigned int &stop, const unsigned int &primo, const unsigned int & ith) noexcept; - std::vector get_shape() const {return {ntimesteps/skip,ntype,3}; } - std::vector get_stride() const {return { static_cast (sizeof (double)*3*ntype),sizeof (double)*3,sizeof (double)};} - void set_starting_center(const std::valarray & s) {if (s.size()!=starting_center.size()) throw std::runtime_error("wrong size of provided starting centers"); starting_center=s;} - void set_zero_threshold(const double & z ){zero=z;} - bool get_reset_error() {bool e=error; error=false; return e;} - void join_data(){} - ~CenterOfMassDiff(); -private: - Traiettoria *t; - using OperazioniSuLista::lista; - using OperazioniSuLista::lunghezza_lista; - unsigned int lista_alloc; - int ntype,maxiter; - std::valarray starting_center; - double zero; - bool error; - -}; - -#endif // CENTEROFMASSDIFF_H diff --git a/lib/include/chargefluxts.h b/lib/include/chargefluxts.h index 12c42f61b..63f5fe50f 100644 --- a/lib/include/chargefluxts.h +++ b/lib/include/chargefluxts.h @@ -12,20 +12,20 @@ #ifndef CHARGEFLUXTS_H #define CHARGEFLUXTS_H -#include "traiettoria.h" +#include "trajectory.h" class ChargeFluxTs { public: - ChargeFluxTs(Traiettoria * t); + ChargeFluxTs(Trajectory * t); ~ChargeFluxTs(); double *J_z(const unsigned int & timestep); - void calcola(unsigned int primo); + void calculate(unsigned int primo); void reset(unsigned int numeroTimestepPerBlocco); private: - Traiettoria * traiettoria; + Trajectory * traiettoria; double *J; - unsigned int timestep_finestra,timestep_corrente; + unsigned int timestep_finestra,current_timestep; bool calcolato; }; diff --git a/lib/include/convolution.h b/lib/include/convolution.h index 700ffc3f4..bee3d58f6 100644 --- a/lib/include/convolution.h +++ b/lib/include/convolution.h @@ -22,8 +22,8 @@ template < typename T >class Convolution **/ Convolution(std::function f, int n, double start, double stop, int center); ~Convolution(); - void calcola(T* in , T *out,int n,int skip); - void calcola(T* in , T *out,int n); + void calculate(T* in , T *out,int n,int skip); + void calculate(T* in , T *out,int n); private: T * fc; unsigned int n_fc,centro; diff --git a/lib/include/correlatorespaziale.h b/lib/include/correlatorespaziale.h index d11ad7791..de39dc735 100644 --- a/lib/include/correlatorespaziale.h +++ b/lib/include/correlatorespaziale.h @@ -12,40 +12,40 @@ #ifndef CORRELATORESPAZIALE_H #define CORRELATORESPAZIALE_H -#include "calcolamultithread.h" +#include "calculatemultithread.h" #include "operazionisulista.h" #include -class Traiettoria; +class Trajectory; -class CorrelatoreSpaziale : public CalcolaMultiThread, public OperazioniSuLista +class CorrelatoreSpaziale : public CalculateMultiThread, public VectorOp { public: - CorrelatoreSpaziale(Traiettoria *t, + CorrelatoreSpaziale(Trajectory *t, std::vector< std::array > k, double sigma2, unsigned int nthreads=0, unsigned int skip=1, bool debug=false ); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b) {return 1;} + unsigned int nExtraTimesteps(unsigned int n_b) {return 1;} void reset(const unsigned int numeroTimestepsPerBlocco); void calc_single_th(const unsigned int &start, const unsigned int &stop, const unsigned int &primo, const unsigned int & ith) noexcept; void s_fac_k(const double k[3], const unsigned int i_t,double * out ) const; int get_sfac_size()const {return size_sfac;} void print(std::ostream & out); - using CalcolaMultiThread::operator=; + using CalculateMultiThread::operator=; ~CorrelatoreSpaziale(); std::vector get_shape() const; std::vector get_stride() const; void join_data(){} private: - using OperazioniSuLista::lista; - using OperazioniSuLista::lunghezza_lista; - Traiettoria *t; + using VectorOp::vdata; + using VectorOp::data_length; + Trajectory *t; double * sfac; std::vector< std::array > klist; double sigma2; diff --git a/lib/include/gofrt.h b/lib/include/gofrt.h index 7ffcb1345..e03f9d41d 100644 --- a/lib/include/gofrt.h +++ b/lib/include/gofrt.h @@ -14,25 +14,25 @@ #define GOFRT_H #include "operazionisulista.h" -#include "calcolamultithread.h" +#include "calculatemultithread.h" namespace Gofrt_Flags { -constexpr int FLAGS = CalcolaMultiThread_Flags::PARALLEL_SPLIT_ATOM | - CalcolaMultiThread_Flags::SERIAL_LOOP_TIME | - CalcolaMultiThread_Flags::SERIAL_LOOP_AVERAGE | - CalcolaMultiThread_Flags::CALL_DEBUG_ROUTINE | - CalcolaMultiThread_Flags::CALL_CALC_INIT; +constexpr int FLAGS = CalculateMultiThread_Flags::PARALLEL_SPLIT_ATOM | + CalculateMultiThread_Flags::SERIAL_LOOP_TIME | + CalculateMultiThread_Flags::SERIAL_LOOP_AVERAGE | + CalculateMultiThread_Flags::CALL_DEBUG_ROUTINE | + CalculateMultiThread_Flags::CALL_CALC_INIT; } template -class Gofrt : public OperazioniSuLista,TFLOAT>, public CalcolaMultiThread, Gofrt_Flags::FLAGS > +class Gofrt : public VectorOp,TFLOAT>, public CalculateMultiThread, Gofrt_Flags::FLAGS > { public: using This = Gofrt; - using CalcolaMultiThread_T = CalcolaMultiThread; - using CalcolaMultiThread_T::FLAGS; - using OperazioniSuLista_T = OperazioniSuLista; + using CalculateMultiThread_T = CalculateMultiThread; + using CalculateMultiThread_T::FLAGS; + using VectorOp_T = VectorOp; Gofrt(T *t, TFLOAT rmin, @@ -45,29 +45,29 @@ class Gofrt : public OperazioniSuLista,TFLOAT>, public CalcolaMu bool debug=false); ~Gofrt(); void reset(const unsigned int numeroTimestepsPerBlocco); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + unsigned int nExtraTimesteps(unsigned int n_b); This & operator =(const This & destra); std::vector get_shape(); std::vector get_stride(); std::string get_columns_description() {return c_descr;} - using OperazioniSuLista_T::azzera; + using VectorOp_T::azzera; void calc_init(int); void calc_single_th(int,int,int,int,int,int); void calc_end(); private: - using OperazioniSuLista_T::lista; - using OperazioniSuLista_T::lunghezza_lista; + using VectorOp_T::vdata; + using VectorOp_T::data_length; TFLOAT * th_data; TFLOAT rmin,rmax,rmax2,rmin2,dr,incr; bool debug; T * traiettoria; unsigned int nbin,lmax; - using CalcolaMultiThread_T::ntimesteps; - using CalcolaMultiThread_T::skip; - using CalcolaMultiThread_T::nthreads; - using CalcolaMultiThread_T::leff; + using CalculateMultiThread_T::ntimesteps; + using CalculateMultiThread_T::skip; + using CalculateMultiThread_T::nthreads; + using CalculateMultiThread_T::leff; int gofr_idx(unsigned int ts, unsigned int itype=0, unsigned int r=0) { unsigned int idx= ts * traiettoria->get_ntypes()*(traiettoria->get_ntypes()+1)*nbin +nbin * itype @@ -76,11 +76,11 @@ class Gofrt : public OperazioniSuLista,TFLOAT>, public CalcolaMu } TFLOAT * gofr(unsigned int ts, unsigned int itype=0, unsigned int r=0){ unsigned int idx= gofr_idx(ts,itype,r); - if (idx >= lunghezza_lista) { + if (idx >= data_length) { std::cerr << "Errore: indice fuori dal range!\n"; abort(); } - return &lista[idx]; + return &vdata[idx]; } std::string c_descr; unsigned int get_itype(unsigned int & type1,unsigned int & type2) const { diff --git a/lib/include/greenkubo2componentionicfluid.h b/lib/include/greenkubo2componentionicfluid.h deleted file mode 100644 index cb53be8c1..000000000 --- a/lib/include/greenkubo2componentionicfluid.h +++ /dev/null @@ -1,50 +0,0 @@ -/** - * - * (c) Riccardo Bertossa, 2019 - * - * Use at your own risk. - * - * If you modified the code, I could be happy if you contribute on github! - * -**/ - -#ifndef GREENKUBO2COMPONENTIONICFLUID_H -#define GREENKUBO2COMPONENTIONICFLUID_H - -#include "readlog.h" - -#include "operazionisulista.h" -#include "mediablocchi.h" -#include - - -class GreenKubo2ComponentIonicFluid : public OperazioniSuLista -{ -public: - GreenKubo2ComponentIonicFluid(ReadLog<> * traiettoria, - std::string log, - double *cariche, - unsigned int skip, - bool dump=false, - unsigned int lunghezza_funzione_max=0, - unsigned int nthreads=0, - unsigned int n_ris=100 - ); - ~GreenKubo2ComponentIonicFluid(); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); - void reset(unsigned int numeroTimestepsPerBlocco); - void calcola(unsigned int primo); - GreenKubo2ComponentIonicFluid & operator =(const GreenKubo2ComponentIonicFluid &); - static const unsigned int narr=14; -private: - bool scrivi_file; - unsigned int idx_je,idx_j0,idx_j1,n_ris; - std::string log; - ReadLog<> *traiettoria; - unsigned int ntimesteps,lmax,leff,nthread,skip; - double carica[2]; - std::array jz(unsigned int ts); - double* je(unsigned int ts); -}; - -#endif // GREENKUBO2COMPONENTIONICFLUID_H diff --git a/lib/include/greenkuboNcomponentionicfluid.h b/lib/include/greenkuboNcomponentionicfluid.h index 9a93af0ff..e796483b5 100644 --- a/lib/include/greenkuboNcomponentionicfluid.h +++ b/lib/include/greenkuboNcomponentionicfluid.h @@ -15,13 +15,12 @@ #include "operazionisulista.h" -#include "mediablocchi.h" #include #include #include template -class GreenKuboNComponentIonicFluid : public OperazioniSuLista,TFLOAT> +class GreenKuboNComponentIonicFluid : public VectorOp,TFLOAT> { public: GreenKuboNComponentIonicFluid(READLOG * traiettoria, @@ -39,9 +38,9 @@ class GreenKuboNComponentIonicFluid : public OperazioniSuLista & operator =(const GreenKuboNComponentIonicFluid &); unsigned int get_narr(); unsigned int get_indexOfKappa(); @@ -52,8 +51,8 @@ class GreenKuboNComponentIonicFluid : public OperazioniSuLista get_stride(); private: - using OperazioniSuLista,TFLOAT>::lista; - using OperazioniSuLista,TFLOAT>::lunghezza_lista; + using VectorOp,TFLOAT>::vdata; + using VectorOp,TFLOAT>::data_length; static bool benchmarked; unsigned int narr,N_corr,start_mean,n_seg,n_seg_start,n_seg_stop; bool scrivi_file,subtract_mean,bench; diff --git a/lib/include/heatc.h b/lib/include/heatc.h deleted file mode 100644 index b3931f715..000000000 --- a/lib/include/heatc.h +++ /dev/null @@ -1,43 +0,0 @@ -/** - * - * (c) Riccardo Bertossa, 2019 - * - * Use at your own risk. - * - * If you modified the code, I could be happy if you contribute on github! - * -**/ - -#ifndef HEATC_H -#define HEATC_H - -#include "operazionisulista.h" -#include - -class Traiettoria; - -class HeatC: public OperazioniSuLista -{ -public: - HeatC(Traiettoria *t, - double sigma, - unsigned int nthreads=0, - unsigned int skip=1); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b) {return 1;} - void reset(const unsigned int numeroTimestepsPerBlocco); - void calcola(unsigned int primo); - HeatC & operator = (const HeatC &); - ~HeatC(); - std::vector get_shape() const; - std::vector get_stride() const; - -private: - double sigma; - unsigned int nthreads; - unsigned int skip; - Traiettoria *t; - unsigned int ntimesteps; - unsigned int tipi_atomi; -}; - -#endif // HEATC_H diff --git a/lib/include/heatfluxts.h b/lib/include/heatfluxts.h index 04069cf2d..50299d20e 100644 --- a/lib/include/heatfluxts.h +++ b/lib/include/heatfluxts.h @@ -14,19 +14,19 @@ #define HEATFLUXTS_H #include -#include "traiettoria.h" +#include "trajectory.h" class HeatFluxTs { public: - HeatFluxTs(std::string filename,Traiettoria * t,unsigned int skip=1); + HeatFluxTs(std::string filename,Trajectory * t,unsigned int skip=1); ~HeatFluxTs(); double * flux(unsigned int ts); double * temp(unsigned int ts); unsigned int get_skip() {return skip;} double get_L(){return L;} private: - Traiettoria * traiettoria; + Trajectory * traiettoria; unsigned int skip; double * heatflux; double * T,L; diff --git a/lib/include/istogrammaatomiraggio.h b/lib/include/istogrammaatomiraggio.h index 19c2340b1..891c1a159 100644 --- a/lib/include/istogrammaatomiraggio.h +++ b/lib/include/istogrammaatomiraggio.h @@ -15,20 +15,20 @@ #include -class Traiettoria; +class Trajectory; class IstogrammaAtomiRaggio { public: - IstogrammaAtomiRaggio(Traiettoria * t,double r,unsigned int skip=1,unsigned int nthreads=0); + IstogrammaAtomiRaggio(Trajectory * t,double r,unsigned int skip=1,unsigned int nthreads=0); void reset(const unsigned int numeroTimestepsPerBlocco); - void calcola(unsigned int); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b) {return 0;} + void calculate(unsigned int); + unsigned int nExtraTimesteps(unsigned int n_b) {return 0;} std::map * get_hist() {return hist;}; private: double r2; unsigned int skip,nthreads,ntypes,natoms,numeroTimestepsBlocco; - Traiettoria * traiettoria; + Trajectory * traiettoria; std::map * hist; }; diff --git a/lib/include/istogrammavelocita.h b/lib/include/istogrammavelocita.h index 4da0edebd..6928bba06 100644 --- a/lib/include/istogrammavelocita.h +++ b/lib/include/istogrammavelocita.h @@ -14,18 +14,18 @@ #define ISTOGRAMMAVELOCITA_H #include "operazionisulista.h" -#include "traiettoria.h" +#include "trajectory.h" -class IstogrammaVelocita : public OperazioniSuLista +class IstogrammaVelocita : public VectorOp { public: - IstogrammaVelocita(Traiettoria *t,unsigned int nbins,double vminmax_); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + IstogrammaVelocita(Trajectory *t,unsigned int nbins,double vminmax_); + unsigned int nExtraTimesteps(unsigned int n_b); void reset(const unsigned int numeroTimestepsPerBlocco); - void calcola(unsigned int primo); + void calculate(unsigned int primo); IstogrammaVelocita & operator =(const IstogrammaVelocita & destra); private: - Traiettoria *traiettoria; + Trajectory *traiettoria; unsigned int ntimestep,bins; double vminmax; diff --git a/lib/include/lammps_struct.h b/lib/include/lammps_struct.h index 21953e092..d6941d7ca 100644 --- a/lib/include/lammps_struct.h +++ b/lib/include/lammps_struct.h @@ -33,7 +33,7 @@ struct Intestazione_timestep { bigint natoms; int triclinic; int condizioni_al_contorno[6]; - double scatola[6]; //xlo,xhi,ylo,yhi,zlo,zhi + double box[6]; //xlo,xhi,ylo,yhi,zlo,zhi int dimensioni_riga_output; int nchunk; char * read(char * file) { @@ -41,7 +41,7 @@ struct Intestazione_timestep { file=read_and_advance(file, &natoms); file=read_and_advance(file,&triclinic); file=read_and_advance(file, condizioni_al_contorno,6); - file=read_and_advance(file,scatola,6); + file=read_and_advance(file,box,6); file=read_and_advance(file,&dimensioni_riga_output); file=read_and_advance(file,&nchunk); return file; @@ -51,7 +51,7 @@ struct Intestazione_timestep { out.write((char*) &natoms, sizeof(bigint)); out.write((char*) &triclinic, sizeof(int)); out.write((char*) condizioni_al_contorno, sizeof(int)*6); - out.write((char*) scatola, sizeof(double)*6); + out.write((char*) box, sizeof(double)*6); out.write((char*)&dimensioni_riga_output, sizeof(int)); out.write((char*)&nchunk, sizeof(int)); } @@ -68,7 +68,7 @@ struct Intestazione_timestep_triclinic { bigint natoms; int triclinic=1; int condizioni_al_contorno[6]; - double scatola[6]; + double box[6]; double xy_xz_yz[3]; int dimensioni_riga_output; int nchunk; @@ -77,7 +77,7 @@ struct Intestazione_timestep_triclinic { file=read_and_advance(file, &natoms); file=read_and_advance(file,&triclinic); file=read_and_advance(file, condizioni_al_contorno,6); - file=read_and_advance(file,scatola,6); + file=read_and_advance(file,box,6); if (triclinic) file=read_and_advance(file,xy_xz_yz,3); file=read_and_advance(file,&dimensioni_riga_output); file=read_and_advance(file,&nchunk); @@ -88,7 +88,7 @@ struct Intestazione_timestep_triclinic { out.write((char*) &natoms, sizeof(bigint)); out.write((char*) &triclinic, sizeof(int)); out.write((char*) condizioni_al_contorno, sizeof(int)*6); - out.write((char*) scatola, sizeof(double)*6); + out.write((char*) box, sizeof(double)*6); if (triclinic) out.write((char*) xy_xz_yz, sizeof(double)*3); out.write((char*)&dimensioni_riga_output, sizeof(int)); out.write((char*)&nchunk, sizeof(int)); @@ -110,7 +110,7 @@ struct Atomo { double id; double tipo; double posizione[3]; - double velocita[3]; + double velocity[3]; static_assert (NDOUBLE_ATOMO==8, "You have to modify class Atomo if you change NDOUBLE_ATOMO" ); static char * read_id_tipo(char * begin, double & id_, double & tipo_) { begin = read_and_advance(begin,&id_); @@ -162,7 +162,7 @@ current_ptr=read_and_advance(current_ptr,what,n); READ_ADV(natoms); READ_ADV(triclinic); READ_ADV_A(condizioni_al_contorno,6); - READ_ADV_A(scatola,6); + READ_ADV_A(box,6); if(triclinic) { READ_ADV_A(xy_xz_yz,3); } @@ -243,7 +243,7 @@ class TimestepManager { select_var(natoms,int) select_var(triclinic,int) select_var(timestep,bigint) - select_var(scatola, double *) + select_var(box, double *) select_var(dimensioni_riga_output,int) double * xy_xz_yz() { diff --git a/lib/include/modivibrazionali.h b/lib/include/modivibrazionali.h index 0beec4fab..8052d332d 100644 --- a/lib/include/modivibrazionali.h +++ b/lib/include/modivibrazionali.h @@ -14,19 +14,20 @@ #define MODIVIBRAZIONALI_H #include "operazionisulista.h" -#include "posizioniequilibrio.h" #include #include +#include "trajectory.h" +#include "posizioniequilibrio.h" -class ModiVibrazionali : public OperazioniSuLista +class ModiVibrazionali : public VectorOp { public: - ModiVibrazionali(Traiettoria * tr, std::string ifcfile, std::string fononefile, unsigned int n_threads, unsigned int timestep_blocco=0); + ModiVibrazionali(Trajectory * tr, std::string ifcfile, std::string fononefile, unsigned int n_threads, unsigned int timestep_blocco=0); void read_force_file(std::string f); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + unsigned int nExtraTimesteps(unsigned int n_b); void azzera(); void reset(unsigned int s); - void calcola(unsigned int primo); + void calculate(unsigned int primo); /* ModiVibrazionali &operator =(const ModiVibrazionali &); @@ -42,7 +43,7 @@ class ModiVibrazionali : public OperazioniSuLista Eigen::Matrix3Xd vettori_onda; std::string fileFononi; unsigned int numero_timesteps,numero_threads,timestepBlocco; - Traiettoria * traiettoria; + Trajectory * traiettoria; PosizioniEquilibrio * posizioni_equilibrio; Eigen::MatrixXd ifc; }; diff --git a/lib/include/mp.h b/lib/include/mp.h index 9a2ba5e87..618ab4ba1 100644 --- a/lib/include/mp.h +++ b/lib/include/mp.h @@ -21,7 +21,7 @@ #include #include -template class OperazioniSuLista; +template class VectorOp; class Mp { @@ -32,12 +32,12 @@ class Mp bool ionode(); std::string outname(std::string s); - template void send_to_root(OperazioniSuLista * l) { - MPI_Send( l->accesso_lista(),l->lunghezza(),MPI_DOUBLE,0,0,MPI_COMM_WORLD); + template void send_to_root(VectorOp * l) { + MPI_Send( l->access_vdata(),l->lunghezza(),MPI_DOUBLE,0,0,MPI_COMM_WORLD); } - template void recv_root(OperazioniSuLista * l,int source){ - MPI_Recv(l->accesso_lista(),l->lunghezza(),MPI_DOUBLE,source,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); + template void recv_root(VectorOp * l,int source){ + MPI_Recv(l->access_vdata(),l->lunghezza(),MPI_DOUBLE,source,0,MPI_COMM_WORLD,MPI_STATUS_IGNORE); } ~Mp(); diff --git a/lib/include/msd.h b/lib/include/msd.h index 81055a309..5dbdd0a37 100644 --- a/lib/include/msd.h +++ b/lib/include/msd.h @@ -14,17 +14,17 @@ #define MSD_H #include "operazionisulista.h" -#include "calcolamultithread.h" +#include "calculatemultithread.h" #include namespace MSD_Flags { -constexpr int FLAGS = CalcolaMultiThread_Flags::PARALLEL_SPLIT_TIME | - CalcolaMultiThread_Flags::CALL_DEBUG_ROUTINE | - CalcolaMultiThread_Flags::CALL_CALC_INIT; +constexpr int FLAGS = CalculateMultiThread_Flags::PARALLEL_SPLIT_TIME | + CalculateMultiThread_Flags::CALL_DEBUG_ROUTINE | + CalculateMultiThread_Flags::CALL_CALC_INIT; } template -class MSD : public OperazioniSuLista >, public CalcolaMultiThread, MSD_Flags::FLAGS > +class MSD : public VectorOp >, public CalculateMultiThread, MSD_Flags::FLAGS > { public: MSD(T *t, @@ -44,12 +44,12 @@ class MSD : public OperazioniSuLista >, public CalcolaMultiThread & operator =(const MSD & destra); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + unsigned int nExtraTimesteps(unsigned int n_b); using This = MSD; - using CMT = CalcolaMultiThread; + using CMT = CalculateMultiThread; private: - using OperazioniSuLista >::lista; - using OperazioniSuLista >::lunghezza_lista; + using VectorOp >::vdata; + using VectorOp >::data_length; T * traiettoria; size_t lmax,f_cm,ntypes; using CMT::nthreads; diff --git a/lib/include/operazionisulista.h b/lib/include/operazionisulista.h index a7f76adca..a952fb6ae 100644 --- a/lib/include/operazionisulista.h +++ b/lib/include/operazionisulista.h @@ -13,185 +13,138 @@ #ifndef OPERAZIONISULISTA_H #define OPERAZIONISULISTA_H -#include "mediablocchi.h" #include "config.h" #include "compiler_types.h" -//#include +#include /** - * Definisce le operazioni, eseguite membro a membro, delle liste. - * Definisce le operazioni, eseguite membro a membro, con uno scalare. + * Simple class to copy and operate on vectors on the heap + * Implements vector-vector elementwise + - * /, and operations with scalars **/ -template class OperazioniSuLista +template class VectorOp { public: - OperazioniSuLista &operator =(const OperazioniSuLista &destra) { -#ifdef DEBUG2 - std::cerr << "chiamato OperazioniSuLista::operator = " __FILE__ ":"<<__LINE__<<"\n"; -#endif - if (lunghezza_lista!=destra.lunghezza_lista) { //rialloca la memoria -#ifdef DEBUG2 - std::cerr << "delete [] "< &operator =(const VectorOp &destra) { + if (data_length!=destra.data_length) { //rialloca la memoria + delete [] vdata; + data_length=destra.data_length; + vdata = new TFLOAT[data_length]; } - if (lunghezza_lista==destra.lunghezza_lista){ - for (unsigned int i=0;i &operator =(OperazioniSuLista &&) = default; + VectorOp &operator =(VectorOp &&) = default; T & operator+= (const T&destra){ - if (destra.lunghezza()== lunghezza_lista) { - for (unsigned int i=0;i(*this); } T & operator-= (const T & destra) { - if (destra.lunghezza()== lunghezza_lista) { - for (unsigned int i=0;i(*this); } T & operator*= (const T & destra) { - if (destra.lunghezza()== lunghezza_lista) { -#ifdef DEBUG2 - std::cerr << "chiamato OperazioniSuLista::operator *= " __FILE__ ":"<<__LINE__<<"\n"; - std::cerr << "lista = "<(*this); } T & operator/= (const T & destra) { - if (destra.lunghezza()== lunghezza_lista) { - for (unsigned int i=0;i(*this); } T & operator+= (const TFLOAT & destra) { - for (unsigned int i=0;i(*this); } T & operator-= (const TFLOAT & destra) { - for (unsigned int i=0;i(*this); } T & operator*= (const TFLOAT & destra) { - for (unsigned int i=0;i(*this); } T & operator/= (const TFLOAT & destra) { - for (unsigned int i=0;i(*this); } -// /**TODO -// * QUESTI NON FUNZIONANO -// * sistemare l'operatore di copia ed risolvere i problemi con la memoria -// **/ -// const T operator+ (const T&) const ; -// const T operator- (const T&) const ; -// const T operator* (const T&) const ; -// const T operator/ (const T&) const ; - -// /**TODO -// * QUESTI NON SONO MAI STATI TESTATI -// * E PROBABILMENTE NON FUNZIONANO -// **/ - -// const T operator+ (const TFLOAT &) const ; -// const T operator- (const TFLOAT&) const ; -// const T operator* (const TFLOAT&) const ; -// const T operator/ (const TFLOAT&) const ; unsigned int lunghezza() const{ - return lunghezza_lista; + return data_length; } TFLOAT elemento(unsigned int i)const{ - if (i (const OperazioniSuLista & copiare) { - lista=0; - lunghezza_lista=0; + VectorOp (const VectorOp & copiare) { + vdata=0; + data_length=0; operator=(copiare); - #ifdef DEBUG - std::cerr << "OperazioniSuLista::OperazioniSuLista(const OperazioniSuLista & copiare) "<::~OperazioniSuLista " __FILE__ ":"<<__LINE__<<"\n"; - std::cerr << "delete [] "< #include #include "eigen_include.h" +#include "trajectory.h" -class PosizioniEquilibrio : public OperazioniSuLista +class PosizioniEquilibrio : public VectorOp { public: - PosizioniEquilibrio(Traiettoria *, unsigned int timesteps_sottoblocco=0); + PosizioniEquilibrio(Trajectory *, unsigned int timesteps_sottoblocco=0); ~PosizioniEquilibrio(); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + unsigned int nExtraTimesteps(unsigned int n_b); void reset(const unsigned int numeroTimestepsPerBlocco); - void calcola(unsigned int primo); + void calculate(unsigned int primo); void fit_nacl(); double d2_reticolo_spostamento_medio(double *min, double *max, double *spostamento); bool zona_brillouin(double * k_test); @@ -60,7 +61,7 @@ class PosizioniEquilibrio : public OperazioniSuLista void reticolo_inizializza(double * base_r, double * base, unsigned int *base_type, const unsigned int & nbase); void coord_reticolo(double *xyz,double *uvw_min,double *uvw_max); bool calcolato; - Traiettoria *traiettoria; + Trajectory *traiettoria; unsigned int lunghezza_media,timestepSottoblocco; }; diff --git a/lib/include/readlog.h b/lib/include/readlog.h index 85a5d80c9..460d2d5c9 100644 --- a/lib/include/readlog.h +++ b/lib/include/readlog.h @@ -13,7 +13,7 @@ #ifndef READLOG_H #define READLOG_H #include -#include "traiettoria.h" +#include "trajectory.h" #include template @@ -21,7 +21,7 @@ class ReadLog { public: ReadLog(std::string filename, - Traiettoria * t=nullptr, + Trajectory * t=nullptr, unsigned int skip=1, unsigned int nthreads=0, unsigned int nbatch=200, @@ -35,13 +35,13 @@ class ReadLog /** * introduco una sintassi strana per calcolare al volo la corrente di carica dal file binario della traiettoria: - * "#traj:JZ N q1 ... qN" --> calcola la corrente dalla traiettoria utilizzando le N cariche per gli atomi q1 ... qN + * "#traj:JZ N q1 ... qN" --> calcola la corrente dalla traiettoria utilizzando le N charge per gli atomi q1 ... qN **/ std::pair get_index_of(std::string header); // int need_binary(std::vector headers); - void calc_currents(Traiettoria * t,unsigned int blocks); + void calc_currents(Trajectory * t,unsigned int blocks); private: - Traiettoria * traiettoria; + Trajectory * traiettoria; std::vector headers; std::vector data; std::vector timesteps; @@ -50,7 +50,7 @@ class ReadLog bool if_only_numbers(std::string str); unsigned int natoms,nthreads; - ///qui vengono memorizzati le stringhe delle correnti da calcolare e le cariche da utilizzare per calcolarle + ///qui vengono memorizzati le stringhe delle correnti da calcolare e le charge da utilizzare per calcolarle std::vector< std::pair > > q_current_type; /// analizza la stringa che definisce la corrente da calcolare std::pair >qs(std::string header); diff --git a/lib/include/spettrovibrazionale.h b/lib/include/spettrovibrazionale.h index cd829fa1a..f03a1d8ac 100644 --- a/lib/include/spettrovibrazionale.h +++ b/lib/include/spettrovibrazionale.h @@ -15,26 +15,24 @@ #include "config.h" -#include "mediablocchi.h" #include "operazionisulista.h" #ifdef HAVEfftw3 #include #else #include #endif -#include "mediablocchi.h" -//#include +#include -class Traiettoria; +class Trajectory; template -class SpettroVibrazionale : public OperazioniSuLista > +class SpettroVibrazionale : public VectorOp > { public: - unsigned int numeroTimestepsOltreFineBlocco(unsigned int n_b); + unsigned int nExtraTimesteps(unsigned int n_b); void reset(const unsigned int numeroTimestepsPerBlocco); - void calcola(unsigned int primo); + void calculate(unsigned int primo); SpettroVibrazionale(T* t,bool dump=false); ~SpettroVibrazionale(); std::vector get_shape() const { return {static_cast (tipi_atomi),static_cast(size/2+1),static_cast(3)} ; } @@ -48,8 +46,8 @@ class SpettroVibrazionale : public OperazioniSuLista > T * traiettoria; unsigned int size; int tipi_atomi; - using OperazioniSuLista >::lista; - using OperazioniSuLista >::lunghezza_lista; + using VectorOp >::vdata; + using VectorOp >::data_length; unsigned int trasformata_size; fftw_complex * trasformata; static fftw_plan fftw3; @@ -60,6 +58,6 @@ class SpettroVibrazionale : public OperazioniSuLista > }; //per fare anche le varie medie a blocchi -//template class MediaBlocchi; +//template class BlockAverage; #endif // SPETTROVIBRAZIONALE_H diff --git a/lib/include/sphericalcorrelations.h b/lib/include/sphericalcorrelations.h index 424e87eea..ec7289372 100644 --- a/lib/include/sphericalcorrelations.h +++ b/lib/include/sphericalcorrelations.h @@ -4,9 +4,10 @@ #include "operazionisulista.h" #include "neighbour.h" #include "sphericalbase.h" +#include template -class SphericalCorrelations : public OperazioniSuLista,TFLOAT>, +class SphericalCorrelations : public VectorOp,TFLOAT>, public SphericalBase { public: @@ -26,10 +27,10 @@ class SphericalCorrelations : public OperazioniSuLista & operator =(const SphericalCorrelations & destra){ - OperazioniSuLista,TFLOAT>::operator = (destra); + VectorOp,TFLOAT>::operator = (destra); return *this; } const std::vector get_shape()const { return {(ssize_t)leff,(ssize_t)ntypes,(ssize_t)ntypes,(ssize_t)nbin,(ssize_t)(l+1)};} @@ -49,7 +50,7 @@ class SphericalCorrelations : public OperazioniSuLista,TFLOAT>::azzera; + using VectorOp,TFLOAT>::azzera; size_t get_single_type_size() const { return SPB::get_single_atom_size(); } @@ -63,8 +64,8 @@ class SphericalCorrelations : public OperazioniSuLista,TFLOAT>::lista; - using OperazioniSuLista,TFLOAT>::lunghezza_lista; + using VectorOp,TFLOAT>::vdata; + using VectorOp,TFLOAT>::data_length; T & t; size_t nbin, tmax,nthreads,skip,leff,ntimesteps,buffer_size; diff --git a/lib/include/steinhardt.h b/lib/include/steinhardt.h index a4161a797..70432a3a4 100644 --- a/lib/include/steinhardt.h +++ b/lib/include/steinhardt.h @@ -2,28 +2,28 @@ #define STEINHARDT_H #include "sphericalbase.h" -#include "calcolamultithread.h" +#include "calculatemultithread.h" #include "operazionisulista.h" namespace Steinhardt_Flags { constexpr int FLAGS = - CalcolaMultiThread_Flags::PARALLEL_SPLIT_AVERAGE | - CalcolaMultiThread_Flags::CALL_INNER_JOIN_DATA | - CalcolaMultiThread_Flags::CALL_CALC_INIT; + CalculateMultiThread_Flags::PARALLEL_SPLIT_AVERAGE | + CalculateMultiThread_Flags::CALL_INNER_JOIN_DATA | + CalculateMultiThread_Flags::CALL_CALC_INIT; } template class Steinhardt : public SphericalBase, - public CalcolaMultiThread< Steinhardt, Steinhardt_Flags::FLAGS >, - public OperazioniSuLista,TFLOAT> + public CalculateMultiThread< Steinhardt, Steinhardt_Flags::FLAGS >, + public VectorOp,TFLOAT> { public: - using LISTA = OperazioniSuLista,TFLOAT>; + using LISTA = VectorOp,TFLOAT>; using SPB = SphericalBase; - using CMT = CalcolaMultiThread, Steinhardt_Flags::FLAGS >; - using CMT::calcola; + using CMT = CalculateMultiThread, Steinhardt_Flags::FLAGS >; + using CMT::calculate; using SPB::calc; using typename SPB::Rminmax_t; using typename SPB::Neighbours_T; @@ -40,7 +40,7 @@ class Steinhardt : public SphericalBase, const bool averaged_order = false ); - unsigned int numeroTimestepsOltreFineBlocco(unsigned int nb) {return 0;} + unsigned int nExtraTimesteps(unsigned int nb) {return 0;} void reset(const unsigned int numeroTimestepsPerBlocco); void calc_init(int primo); void calc_single_th(int,//average index, begin @@ -110,8 +110,8 @@ class Steinhardt : public SphericalBase, const size_t nbin; std::string c_descr; T & t; - using LISTA::lunghezza_lista; - using LISTA::lista; + using LISTA::data_length; + using LISTA::vdata; const std::vector steinhardt_l_histogram; size_t nbin_steinhardt,steinhardt_histogram_size; TFLOAT * threadResults; diff --git a/lib/include/testtraiettoria.h b/lib/include/testtraiettoria.h index 26fc7ee3b..e7dee0d0a 100644 --- a/lib/include/testtraiettoria.h +++ b/lib/include/testtraiettoria.h @@ -13,11 +13,11 @@ #ifndef TESTTRAIETTORIA_H #define TESTTRAIETTORIA_H -#include "traiettoria.h" +#include "trajectory.h" -class Traiettoria; +class Trajectory; -class TestTraiettoria : public Traiettoria +class TestTraiettoria : public Trajectory { public: TestTraiettoria(std::string filename); diff --git a/lib/include/traiettoria_numpy.h b/lib/include/traiettoria_numpy.h deleted file mode 100644 index 672612e75..000000000 --- a/lib/include/traiettoria_numpy.h +++ /dev/null @@ -1,38 +0,0 @@ -#ifndef TRAIETTORIA_NUMPY_H -#define TRAIETTORIA_NUMPY_H - -#include "traiettoriabase.h" -#include "pybind11/pybind11.h" - -class Traiettoria_numpy : public TraiettoriaBase -{ -public: - Traiettoria_numpy(pybind11::buffer buffer_pos, - pybind11::buffer buffer_vel, - pybind11::buffer buffer_types, - pybind11::buffer buffer_box, ///data of the boxes, in the format specified by matrix_box - TraiettoriaBase::BoxFormat matrix_box=TraiettoriaBase::BoxFormat::Cell_vectors, - bool pbc_wrap=false, - bool save_rotation_matrix=false); - ~Traiettoria_numpy(); - template - double * posizioni (const int & timestep, const int & atomo) {return buffer_posizioni+natoms*3*timestep+atomo*3;} - template - double * velocita (const int & timestep, const int & atomo){return buffer_velocita+natoms*3*timestep+atomo*3;} - template - double * scatola (const int & timestep){return buffer_scatola+timestep*buffer_scatola_stride;} - template - double * posizioni_cm(const int & timestep, const int & tipo){return buffer_posizioni_cm+timestep*ntypes*3 + tipo*3;} - template - double * velocita_cm(const int & timestep, const int & tipo){return buffer_velocita_cm+timestep*ntypes*3 + tipo*3;} - double *scatola_last(){return buffer_scatola + (n_timesteps-1)*buffer_scatola_stride; } - double * get_rotation_matrix(size_t t){return rotation_matrix+9*t;} - using TraiettoriaBase::BoxFormat; -private: - pybind11::buffer buffer_pos,buffer_vel,buffer_types,buffer_box; - double * rotation_matrix; - bool posizioni_allocated, velocities_allocated, box_allocated; - void calc_cm_pos_vel(double * a, double *&cm); -}; - -#endif // TRAIETTORIA_NUMPY_H diff --git a/lib/include/traiettoria.h b/lib/include/trajectory.h similarity index 67% rename from lib/include/traiettoria.h rename to lib/include/trajectory.h index f609a001c..adaeba5d4 100644 --- a/lib/include/traiettoria.h +++ b/lib/include/trajectory.h @@ -20,25 +20,25 @@ #include #include -#include "traiettoriabase.h" +#include "basetrajectory.h" #include "lammps_struct.h" /** - * Questa classe garantisce un accesso veloce ai timestep caricati con imposta_inizio_accesso + * Questa classe garantisce un accesso veloce ai timestep caricati con set_access_at * , che carica tutti i timestepa partendo da quello specificato in un numero pari a quello richiesto - * con la funzione imposta_dimensione_finestra_accesso in precedenza. Caricare nuovi pezzi costa nuove + * con la funzione set_data_access_block_size in precedenza. Caricare nuovi pezzi costa nuove * allocazioni (e deallocazioni) di memoria e tempo cpu. Il file viene letto tramite la chiamata di sistema mmap. * Nel programma vengono allocati con new [] solo degli array dove vengono immagazzinati i dati * della finestra (quindi consumando meno memoria che nel caricamento del file in un unico colpo) **/ -class Traiettoria : public TraiettoriaBase +class Trajectory : public BaseTrajectory { public: - Traiettoria(std::string filename); - ~Traiettoria(); + Trajectory(std::string filename); + ~Trajectory(); template double * get_array(double * base, const size_t & timestep, const size_t & atomo, const size_t & stride1, const size_t &stride2) { @@ -51,17 +51,17 @@ class Traiettoria : public TraiettoriaBase } } - size_t t=timestep-timestep_corrente; // note that this can be garbage if timestep=timestep_corrente) { // vuol dire che ho già caricato i dati + if (dati_caricati && t < loaded_timesteps && timestep>=current_timestep) { // vuol dire che ho già caricato i dati return base+t*stride1+atomo*stride2; } else { // non ho caricato i dati, li carico prima (questo potrebbe essere inefficiente se dopo devo satare di nuovo indietro! #ifdef DEBUG std::cerr << "Warning: loading timesteps starting from "<=timestep_corrente && t< loaded_timesteps){ + if(set_access_at(timestep)) { + t=timestep-current_timestep; + if (timestep>=current_timestep && t< loaded_timesteps){ return base + t*stride1+atomo*stride2; }else{ throw std::runtime_error("requested timestep is out of range"); @@ -71,42 +71,42 @@ class Traiettoria : public TraiettoriaBase } } } else { - return base + (timestep-timestep_corrente)*stride1+atomo*stride2; + return base + (timestep-current_timestep)*stride1+atomo*stride2; } } template - double * posizioni (const size_t ×tep, const size_t &atomo){ - return get_array(buffer_posizioni,timestep,atomo,3*natoms,3); + double * positions (const size_t ×tep, const size_t &atomo){ + return get_array(buffer_positions,timestep,atomo,3*natoms,3); } template - double * velocita(const size_t ×tep, const size_t &atomo) { - return get_array(buffer_velocita,timestep,atomo,3*natoms,3); + double * velocity(const size_t ×tep, const size_t &atomo) { + return get_array(buffer_velocity,timestep,atomo,3*natoms,3); } template - double * scatola(const size_t ×tep) { - return get_array(buffer_scatola,timestep,0,buffer_scatola_stride,0); + double * box(const size_t ×tep) { + return get_array(buffer_boxes,timestep,0,buffer_boxes_stride,0); } template - double * posizioni_cm(const size_t ×tep, const size_t &tipo){ - return get_array(buffer_posizioni_cm,timestep,tipo,3*ntypes,3); + double * positions_cm(const size_t ×tep, const size_t &tipo){ + return get_array(buffer_positions_cm,timestep,tipo,3*ntypes,3); } template - double * velocita_cm(const size_t ×tep, const size_t &tipo){ - return get_array(buffer_velocita_cm,timestep,tipo,3*ntypes,3); + double * velocity_cm(const size_t ×tep, const size_t &tipo){ + return get_array(buffer_velocity_cm,timestep,tipo,3*ntypes,3); } - double *scatola_last(){ + double *box_last(){ if (loaded_timesteps>0) { - return &buffer_scatola[ buffer_scatola_stride*(loaded_timesteps-1)]; + return &buffer_boxes[ buffer_boxes_stride*(loaded_timesteps-1)]; } else { throw std::runtime_error("No data is loaded!\n"); } } - using TraiettoriaBase::Errori; - Traiettoria::Errori imposta_dimensione_finestra_accesso(const size_t & timesteps); - Traiettoria::Errori imposta_inizio_accesso(const size_t & timesteps); + using BaseTrajectory::Errori; + Trajectory::Errori set_data_access_block_size(const size_t & timesteps); + Trajectory::Errori set_access_at(const size_t & timesteps); int64_t get_timestep_lammps(size_t timestep); void index_all(); int * get_lammps_id(); diff --git a/lib/include/trajectory_numpy.h b/lib/include/trajectory_numpy.h new file mode 100644 index 000000000..99bfd0182 --- /dev/null +++ b/lib/include/trajectory_numpy.h @@ -0,0 +1,38 @@ +#ifndef TRAIETTORIA_NUMPY_H +#define TRAIETTORIA_NUMPY_H + +#include "basetrajectory.h" +#include "pybind11/pybind11.h" + +class Trajectory_numpy : public BaseTrajectory +{ +public: + Trajectory_numpy(pybind11::buffer buffer_pos, + pybind11::buffer buffer_vel, + pybind11::buffer buffer_types, + pybind11::buffer buffer_box, ///data of the boxes, in the format specified by matrix_box + BaseTrajectory::BoxFormat matrix_box=BaseTrajectory::BoxFormat::Cell_vectors, + bool pbc_wrap=false, + bool save_rotation_matrix=false); + ~Trajectory_numpy(); + template + double * positions (const int & timestep, const int & atomo) {return buffer_positions+natoms*3*timestep+atomo*3;} + template + double * velocity (const int & timestep, const int & atomo){return buffer_velocity+natoms*3*timestep+atomo*3;} + template + double * box (const int & timestep){return buffer_boxes+timestep*buffer_boxes_stride;} + template + double * positions_cm(const int & timestep, const int & tipo){return buffer_positions_cm+timestep*ntypes*3 + tipo*3;} + template + double * velocity_cm(const int & timestep, const int & tipo){return buffer_velocity_cm+timestep*ntypes*3 + tipo*3;} + double *box_last(){return buffer_boxes + (n_timesteps-1)*buffer_boxes_stride; } + double * get_rotation_matrix(size_t t){return rotation_matrix+9*t;} + using BaseTrajectory::BoxFormat; +private: + pybind11::buffer buffer_pos,buffer_vel,buffer_types,buffer_box; + double * rotation_matrix; + bool posizioni_allocated, velocities_allocated, box_allocated; + void calc_cm_pos_vel(double * a, double *&cm); +}; + +#endif // TRAIETTORIA_NUMPY_H diff --git a/lib/src/atomicdensity.cpp b/lib/src/atomicdensity.cpp index 7b946fb4a..70d5f22f4 100644 --- a/lib/src/atomicdensity.cpp +++ b/lib/src/atomicdensity.cpp @@ -3,17 +3,17 @@ template AtomicDensity::AtomicDensity(T *t, std::array nbin, unsigned int nthreads, unsigned int skip) : - CalcolaMultiThread {nthreads, skip}, nbin{nbin},t{t}, ntypes{t->get_ntypes()} + CalculateMultiThread {nthreads, skip}, nbin{nbin},t{t}, ntypes{t->get_ntypes()} { - lunghezza_lista=nbin[0]*nbin[1]*nbin[2]*ntypes; - lista=new Hist [lunghezza_lista]; - hist=new Hist [ lunghezza_lista*(nthreads-1)]; + data_length=nbin[0]*nbin[1]*nbin[2]*ntypes; + vdata=new Hist [data_length]; + hist=new Hist [ data_length*(nthreads-1)]; azzera(); } template AtomicDensity::~AtomicDensity() { - delete [] hist; // lista deleted in OperazioniSuLista + delete [] hist; // vdata deleted in VectorOp } template @@ -22,14 +22,14 @@ void AtomicDensity::calc_single_th(const unsigned int &start, const uns throw std::runtime_error("trajectory is too short for this kind of calculation. Select a different starting timestep or lower the size of the block"); } Hist * hist_=0; - if (ith>0) hist_= hist+lunghezza_lista*(ith-1); - else hist_= lista; + if (ith>0) hist_= hist+data_length*(ith-1); + else hist_= vdata; unsigned int natoms=t->get_natoms(); for (unsigned int i=start;iget_type(iatom); - auto * pos = t->posizioni(i+primo,iatom); - auto * l = t->scatola(i+primo); + auto * pos = t->positions(i+primo,iatom); + auto * l = t->box(i+primo); bool in_range=false; auto ih=idx_(pos,l,in_range,itype); if (in_range){ @@ -43,19 +43,19 @@ void AtomicDensity::calc_single_th(const unsigned int &start, const uns template void AtomicDensity::join_data() { for (unsigned int ith=0;ith -unsigned int AtomicDensity::numeroTimestepsOltreFineBlocco(unsigned int n_b) {return 0;} +unsigned int AtomicDensity::nExtraTimesteps(unsigned int n_b) {return 0;} template void AtomicDensity::reset(const unsigned int numeroTimestepsPerBlocco) { ntimesteps=numeroTimestepsPerBlocco; azzera(); - for (int i=0;i @@ -68,13 +68,13 @@ std::vector AtomicDensity::get_stride() const { } #ifdef BUILD_MMAP -#include "traiettoria.h" -template class AtomicDensity; +#include "trajectory.h" +template class AtomicDensity; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class AtomicDensity; +#include "trajectory_numpy.h" +template class AtomicDensity; #endif diff --git a/lib/src/traiettoriabase.cpp b/lib/src/basetrajectory.cpp similarity index 71% rename from lib/src/traiettoriabase.cpp rename to lib/src/basetrajectory.cpp index de96b820f..d829b0dce 100644 --- a/lib/src/traiettoriabase.cpp +++ b/lib/src/basetrajectory.cpp @@ -1,10 +1,10 @@ #include "testtraiettoria.h" -#include "traiettoria.h" -#include "traiettoria_numpy.h" +#include "trajectory.h" +#include "trajectory_numpy.h" #include "config.h" template -void TraiettoriaBase::dump_lammps_bin_traj(const std::string &fname, int start_ts, int stop_ts){ +void BaseTrajectory::dump_lammps_bin_traj(const std::string &fname, int start_ts, int stop_ts){ if (start_ts<0 || start_ts>=n_timesteps){ throw std::runtime_error("You must provide a starting timestep between 0 and the number of timesteps!"); } @@ -15,13 +15,13 @@ void TraiettoriaBase::dump_lammps_bin_traj(const std::string &fname, int star Intestazione_timestep_triclinic head; head.natoms=natoms; for (unsigned int i=0;i<6;++i) - head.scatola[i]=scatola(t)[i]; - internal_to_lammps(head.scatola); + head.box[i]=box(t)[i]; + internal_to_lammps(head.box); head.timestep=t; head.triclinic=triclinic; if (triclinic) { for (unsigned int i=0;i<3;++i) - head.xy_xz_yz[i]=scatola(t)[6+i]; + head.xy_xz_yz[i]=box(t)[6+i]; } head.condizioni_al_contorno[0]=0; head.condizioni_al_contorno[1]=0; @@ -41,8 +41,8 @@ void TraiettoriaBase::dump_lammps_bin_traj(const std::string &fname, int star data[0]=iatom; data[1]=get_type(iatom); for (int i=0;i<3;++i){ - data[2+i]=posizioni(t,iatom)[i]; - data[5+i]=velocita(t,iatom)[i]; + data[2+i]=positions(t,iatom)[i]; + data[5+i]=velocity(t,iatom)[i]; } static_assert (NDOUBLE_ATOMO==8, "You have to change the file writing (what do I have to write?) this if you change NDOUBLE_ATOMO" ); out.write((char*) data,sizeof(double)*NDOUBLE_ATOMO); @@ -51,26 +51,26 @@ void TraiettoriaBase::dump_lammps_bin_traj(const std::string &fname, int star } template -size_t TraiettoriaBase::get_ntypes (){ +size_t BaseTrajectory::get_ntypes (){ if (ntypes==0) { types.clear(); - min_type=buffer_tipi[0]; - max_type=buffer_tipi[0]; + min_type=buffer_type[0]; + max_type=buffer_type[0]; bool *duplicati = new bool[natoms]; for (size_t i=0;imax_type) - max_type=buffer_tipi[i]; - if (buffer_tipi[i]max_type) + max_type=buffer_type[i]; + if (buffer_type[i]::get_ntypes (){ type_map[types[i]]=i; } delete [] duplicati; - masse = new double [ntypes]; - cariche = new double [ntypes]; + mass = new double [ntypes]; + charge = new double [ntypes]; for (size_t i=0;i; +template class BaseTrajectory; #endif #ifdef BUILD_MMAP -template class TraiettoriaBase; +template class BaseTrajectory; #endif diff --git a/lib/src/centerdiff.cpp b/lib/src/centerdiff.cpp deleted file mode 100644 index 130d41f92..000000000 --- a/lib/src/centerdiff.cpp +++ /dev/null @@ -1,84 +0,0 @@ -#include "centerdiff.h" -#include -#include - -CenterDiff::CenterDiff(Traiettoria *t, unsigned int nthreads, unsigned int skip, unsigned int nit, bool sum_first_two_and_ignore_vz,bool sum_1) : - CalcolaMultiThread (nthreads,skip), t{t}, nit{nit},lista_alloc{0},starting_center{},sum_first_two_and_ignore_vz{sum_first_two_and_ignore_vz},sum_1{sum_1} -{ - if (nit==0) nit=1; -} - -void CenterDiff::reset(const unsigned int numeroTimestepsPerBlocco) { - lunghezza_lista=numeroTimestepsPerBlocco*3*3*nit/skip; - ntimesteps=numeroTimestepsPerBlocco; - if (lunghezza_lista> lista_alloc){ - delete [] lista; - lista = new double[lunghezza_lista]; - lista_alloc=lunghezza_lista; - } -} - - - -void CenterDiff::calc_single_th(const unsigned int &start, const unsigned int &stop, const unsigned int &primo, const unsigned int &ith) noexcept { - unsigned int ilista=(start-primo)/skip; - std::array dx,cn; - std::array sums; - dx=starting_center; - for (unsigned int itimestep=start;itimestepscatola(itimestep); - for (unsigned int i=0;iget_natoms();++iatom) { - double * coord = t->posizioni(itimestep,iatom); - double * vel = t->velocita(itimestep,iatom); - for (unsigned int idim1=0;idim1<3;idim1++){ - double l=box[idim1*2+1]-box[idim1*2+0]; - if (sum_first_two_and_ignore_vz && idim1==2){ - if (sum_1){ - sums[idim1]+=1.0; - } else { - sums[idim1]+=vel[0]+vel[1]; - } - } else { - sums[idim1]+=vel[idim1]; - } - if (!sum_first_two_and_ignore_vz){ - for (unsigned int idim2=0;idim2<3;idim2++) { - double new_coord=coord[idim1]-dx[3*idim2+idim1]; //shift of the origin - double wrapped_coord=new_coord-std::round(new_coord/l)*l; // center of the cell is in zero - cn[idim2*3+idim1]+=wrapped_coord*vel[idim2]; - } - } else { - for (unsigned int idim2=0;idim2<2;idim2++) { - double new_coord=coord[idim1]-dx[3*idim2+idim1]; //shift of the origin - double wrapped_coord=new_coord-std::round(new_coord/l)*l; // center of the cell is in zero - cn[idim2*3+idim1]+=wrapped_coord*vel[idim2]; - } - unsigned int idim2=2; - double new_coord=coord[idim1]-dx[3*idim2+idim1]; //shift of the origin - double wrapped_coord=new_coord-std::round(new_coord/l)*l; // center of the cell is in zero - if (sum_1){ - cn[idim2*3+idim1]+=wrapped_coord; - } else { - cn[idim2*3+idim1]+=wrapped_coord*(vel[0]+vel[1]); - } - - } - } - } - for (unsigned int idim1=0;idim1<3;++idim1){ - for (unsigned int idim2=0;idim2<3;++idim2){ - dx[idim2*3+idim1]=dx[idim2*3+idim1]+cn[idim2*3+idim1]/sums[idim2]; - lista[ilista*3*3*nit+i*3*3+idim2*3+idim1]=dx[idim2*3+idim1]; - } - } - } - ++ilista; - } -} - -CenterDiff::~CenterDiff() { - -} diff --git a/lib/src/centerofmassdiff.cpp b/lib/src/centerofmassdiff.cpp deleted file mode 100644 index 32c5ac3cd..000000000 --- a/lib/src/centerofmassdiff.cpp +++ /dev/null @@ -1,69 +0,0 @@ -#include "centerofmassdiff.h" - -CenterOfMassDiff::CenterOfMassDiff(Traiettoria *t, unsigned int nthreads, unsigned int skip) : - CalcolaMultiThread (nthreads,skip), t{t},lista_alloc{0},starting_center{},ntype{0},zero{1.0e-10},maxiter{100},error{false} -{ - -} - -void CenterOfMassDiff::reset(const unsigned int numeroTimestepsPerBlocco) { - ntype=t->get_ntypes(); - ntimesteps=numeroTimestepsPerBlocco; - lunghezza_lista=ntimesteps*3*ntype/skip; - starting_center.resize(3*ntype); - if (lunghezza_lista>lista_alloc) { - delete [] lista; - lista = new double [lunghezza_lista]; - lista_alloc=lunghezza_lista; - } -} - -void CenterOfMassDiff::calc_single_th(const unsigned int &start, const unsigned int &stop, const unsigned int &primo, const unsigned int &ith) noexcept { - unsigned int ilista=(start-primo)/skip; - std::valarray dx,cn; - std::valarray sums; - dx=starting_center; - cn.resize((starting_center.size())); - sums.resize(ntype); - for (unsigned int itimestep=start;itimestepscatola(itimestep); - double cns=0.0; - int iter=0; - do { - sums=0.0; - cn=0.0; - for (unsigned int iatom=0;iatomget_natoms();++iatom) { - double * coord = t->posizioni(itimestep,iatom); - unsigned int itype=t->get_type(iatom); - sums[itype]++; - for (unsigned int idim1=0;idim1<3;idim1++){ - double l=box[idim1*2+1]-box[idim1*2+0]; - double new_coord=coord[idim1]-dx[3*itype+idim1]; //shift of the origin, according to the old center - double wrapped_coord=new_coord-std::round(new_coord/l)*l; // center of the cell is in zero - cn[itype*3+idim1]+=wrapped_coord; - } - } - cns=0.0; - for (unsigned int itype=0;itypezero && ++iter=maxiter) { - error=true; - } - for (unsigned int itype=0;itype0){ @@ -76,21 +76,21 @@ void ChargeFluxTs::calcola(unsigned int timestep) { for (unsigned int iatom=0;iatomget_natoms();iatom++) { for (unsigned int i=0;i<3;i++){ - J[istep*3+i]+=traiettoria->velocita(t,iatom)[i] + J[istep*3+i]+=traiettoria->velocity(t,iatom)[i] *traiettoria->get_charge(traiettoria->get_type(iatom)); } } for (unsigned int i=0;i<3;i++) J[istep*3+i]/=double(traiettoria->get_natoms()); } - timestep_corrente=timestep; + current_timestep=timestep; calcolato=true; } double * ChargeFluxTs::J_z(const unsigned int & timestep) { - if (timestep >= timestep_corrente && timestep = current_timestep && timestep > head.natoms >> head.scatola[0] >> head.scatola[1] >> head.scatola[2] >> head.scatola[3] >> head.scatola[4] >> head.scatola[5]; + in >> head.natoms >> head.box[0] >> head.box[1] >> head.box[2] >> head.box[3] >> head.box[4] >> head.box[5]; head.timestep=itim; if (!in.good()){ if (in.eof()) { @@ -116,20 +116,20 @@ ConvertiBinario::ConvertiBinario(const std::string filein, const std::string fil std::cerr << "Fine del file raggiunta dopo "< Convolution::Convolution(std::function void Convolution::calcola(T *in, T*out, int n, int skip) { +template void Convolution::calculate(T *in, T*out, int n, int skip) { for (unsigned int in_index=0;in_index void Convolution::calcola(T *in, T*out, int n, int skip } } -template void Convolution::calcola(T *in, T*out, int n) { +template void Convolution::calculate(T *in, T*out, int n) { for (unsigned int in_index=0;in_index #include "correlatorespaziale.h" -#include "traiettoria.h" +#include "trajectory.h" #include "operazionisulista.h" #include "config.h" #define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628 -CorrelatoreSpaziale::CorrelatoreSpaziale(Traiettoria *t, +CorrelatoreSpaziale::CorrelatoreSpaziale(Trajectory *t, std::vector< std::array > ks, double sigma2, unsigned int nthreads, unsigned int skip, bool debug) : - t{t},sfac(nullptr), sigma2(sigma2),debug(debug),tipi_atomi(0),ntimesteps(0), CalcolaMultiThread (nthreads,skip) + t{t},sfac(nullptr), sigma2(sigma2),debug(debug),tipi_atomi(0),ntimesteps(0), CalculateMultiThread (nthreads,skip) { nk=ks.size(); klist=ks; @@ -42,14 +42,14 @@ void CorrelatoreSpaziale::reset(const unsigned int numeroTimestepsPerBlocco) { return; } delete [] sfac; - delete [] lista; + delete [] vdata; tipi_atomi=t->get_ntypes(); ntimesteps=numeroTimestepsPerBlocco; size_k=3*tipi_atomi*2;// (complex number) size_sfac=size_k*nk; - lunghezza_lista=ntimesteps*size_sfac/skip; + data_length=ntimesteps*size_sfac/skip; sfac = new double[size_sfac*nthreads]; - lista = new double[lunghezza_lista]; + vdata = new double[data_length]; } @@ -70,8 +70,8 @@ void CorrelatoreSpaziale::s_fac_k(const double k[3], double *out1 ///size 3*type {(x,y,z)_1, ..., (x,y,z)_ntype} ) const { for (unsigned int i_at=0;i_atget_natoms();i_at++){ - const double * pos = t->posizioni(i_t,i_at); - const double * vel = t->velocita(i_t,i_at); + const double * pos = t->positions(i_t,i_at); + const double * vel = t->velocity(i_t,i_at); unsigned int type=t->get_type(i_at); double arg=k[0]*pos[0]+k[1]*pos[1]+k[2]*pos[2]; double s=sin(arg),c=cos(arg); @@ -87,7 +87,7 @@ void CorrelatoreSpaziale::calc_single_th(const unsigned int &start, const unsign double * sfact=&sfac[size_sfac*ith]; int ilista=(start-primo)/skip; for (unsigned int itimestep=start;itimestep #include "config.h" -#ifdef USE_MPI -#include "mp.h" -#endif - template Gofrt::Gofrt(T *t, TFLOAT rmin, TFLOAT rmax, unsigned int nbin, unsigned int tmax, unsigned int nthreads, unsigned int skip,unsigned int every, bool debug) : - CalcolaMultiThread_T {nthreads, skip, t->get_natoms(), every}, + CalculateMultiThread_T {nthreads, skip, t->get_natoms(), every}, traiettoria(t),rmin(rmin),rmax(rmax),nbin(nbin), lmax(tmax), debug(debug) { @@ -37,7 +33,7 @@ template Gofrt::~Gofrt() { } -template unsigned int Gofrt::numeroTimestepsOltreFineBlocco(unsigned int n_b){ +template unsigned int Gofrt::nExtraTimesteps(unsigned int n_b){ return (traiettoria->get_ntimesteps()/(n_b+1)+1 < lmax || lmax==0)? traiettoria->get_ntimesteps()/(n_b+1)+1 : lmax; } @@ -59,10 +55,10 @@ template void Gofrt::reset(const unsigned int leff =(numeroTimestepsPerBloccoget_ntypes()*(traiettoria->get_ntypes()+1)*nbin; + data_length=leff*traiettoria->get_ntypes()*(traiettoria->get_ntypes()+1)*nbin; - delete [] lista; - lista=new TFLOAT [lunghezza_lista]; + delete [] vdata; + vdata=new TFLOAT [data_length]; } template std::vector Gofrt::get_shape(){ @@ -86,9 +82,9 @@ void Gofrt::calc_init(int primo) { } //init - th_data = new TFLOAT[lunghezza_lista*(nthreads-1)]; + th_data = new TFLOAT[data_length*(nthreads-1)]; azzera(); - for (unsigned int i=0;i0) incr=1.0/int(ntimesteps/skip); @@ -97,9 +93,9 @@ void Gofrt::calc_init(int primo) { template void Gofrt::calc_single_th(int t, int imedia, int atom_start, int atom_stop,int primo, int ith) { - TFLOAT * th_data_ = th_data + (ith-1)*lunghezza_lista; + TFLOAT * th_data_ = th_data + (ith-1)*data_length; if (ith==0) { - th_data_ = lista; + th_data_ = vdata; } for (unsigned int iatom=atom_start;iatomget_natoms();jatom++) { @@ -129,8 +125,8 @@ void Gofrt::calc_single_th(int t, int imedia, int atom_start, int atom template void Gofrt::calc_end() { for (int ith=0;ith::calc_end() { } template Gofrt & Gofrt::operator =(const Gofrt &destra) { - OperazioniSuLista,TFLOAT >::operator = (destra); + VectorOp,TFLOAT >::operator = (destra); return *this; } #ifdef BUILD_MMAP -template class Gofrt; +#include "trajectory.h" +template class Gofrt; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class Gofrt; +#include "trajectory_numpy.h" +template class Gofrt; #endif diff --git a/lib/src/greenkubo2componentionicfluid.cpp b/lib/src/greenkubo2componentionicfluid.cpp deleted file mode 100644 index 350835dae..000000000 --- a/lib/src/greenkubo2componentionicfluid.cpp +++ /dev/null @@ -1,305 +0,0 @@ -/** - * - * (c) Riccardo Bertossa, 2019 - * - * Use at your own risk. - * - * If you modified the code, I could be happy if you contribute on github! - * -**/ - -#include "greenkubo2componentionicfluid.h" -#include "chargefluxts.h" -#include "heatfluxts.h" -#include -#include -#include -#include -#include - -//const unsigned int GreenKubo2ComponentIonicFluid::narr=14; - -GreenKubo2ComponentIonicFluid::GreenKubo2ComponentIonicFluid(ReadLog<> *traiettoria, std::string log, double * cariche, unsigned int skip, bool dump, unsigned int lunghezza_funzione_max, unsigned int nthreads,unsigned int n_ris) : OperazioniSuLista(), - traiettoria (traiettoria), log(log), ntimesteps(0),skip(skip), scrivi_file(dump),lmax(lunghezza_funzione_max),nthread(nthreads),n_ris(n_ris) -{ - carica[0]=cariche[0]; - carica[1]=cariche[1]; - - std::pair res; - res=traiettoria->get_index_of("c_flusso[1]"); - if(res.second) - idx_je=res.first; - else - {std::cerr << "Errore: header di una corrente non trovato!";assert(0);abort();} - res=traiettoria->get_index_of("c_vcm[1][1]"); - if(res.second) - idx_j0=res.first; - else - {std::cerr << "Errore: header di una corrente non trovato!";assert(0);abort();} - res=traiettoria->get_index_of("c_vcm[2][1]"); - if(res.second) - idx_j1=res.first; - else - {std::cerr << "Errore: header di una corrente non trovato!";assert(0);abort();} - -} - -GreenKubo2ComponentIonicFluid::~GreenKubo2ComponentIonicFluid(){ -#ifdef DEBUG2 -// std::cerr << "Called delete GreenKubo2ComponentIonicFluid, je="<::operator =( destra); - return *this; -} - -unsigned int GreenKubo2ComponentIonicFluid::numeroTimestepsOltreFineBlocco(unsigned int n_b) { - return (traiettoria->n_timestep()/(n_b+1)+1 < lmax || lmax==0)? traiettoria->n_timestep()/(n_b+1)+1 : lmax; -} - -void GreenKubo2ComponentIonicFluid::reset(unsigned int numeroTimestepsPerBlocco) { - leff=(numeroTimestepsPerBlocco GreenKubo2ComponentIonicFluid::jz(unsigned int ts){ - double *j0=&traiettoria->line(ts)[idx_j0]; - double *j1=&traiettoria->line(ts)[idx_j1]; - return std::array {{ - j0[0]*carica[0]+j1[0]*carica[1], - j0[1]*carica[0]+j1[1]*carica[1], - j0[2]*carica[0]+j1[2]*carica[1], - }}; -} - -double* GreenKubo2ComponentIonicFluid::je(unsigned int ts){ - return&traiettoria->line(ts)[idx_je]; -} - -void GreenKubo2ComponentIonicFluid::calcola(unsigned int primo) { - - - unsigned int allinea=0;//primo%skip; - - double intzz=0.0; - double intee=0.0; - double intez=0.0; - double intze=0.0; - double jeeo=0.0,jzzo=0.0,jezo=0.0,jzeo=0.0; - double int_ein_ee=0.0,int_ein_ze=0.0,int_ein_ez=0.0,int_ein_zz=0.0; - - - if (nthread<1){ - std::cerr << "Sorry, not implemented. Please specify '-N #thread number' option.\n"; - abort(); - for (unsigned int itimestep=0;itimestep 1 - - /*dividi il lavoro in gruppi - itimestep ---> [0,leff[ - npassi ---> leff - - ciascun gruppo avrà npassith=npassi/nthread passi - l'ultimo deve finire alla fine - - gli integrali li calcolo alla fine, è uguale. - Nel caso di 1 thread tengo l'implementazione precedente - (anche per poter confrontare i risultati) - */ - - unsigned int npassith=leff/nthread; - std::vector threads; - for (unsigned int ith=0;ith1) { - for (unsigned int itimestep=npassith;itimestep GreenKuboNComponentIoni unsigned int n_seg, bool do_bench, unsigned int n_seg_start, - unsigned int n_seg_stop) : OperazioniSuLista,TFLOAT>(), + unsigned int n_seg_stop) : VectorOp,TFLOAT>(), traiettoria (traiettoria), log(log), ntimesteps(0),skip(skip), scrivi_file(dump), lmax(lunghezza_funzione_max),nthread(nthreads),subtract_mean(subtract_mean), start_mean(start_mean),n_seg(n_seg),bench(false), @@ -111,22 +111,22 @@ template GreenKuboNComponentIoni #ifdef DEBUG2 std::cerr << "Chiamato GreenKuboNComponentIonicFluid::operator =\n"; #endif - OperazioniSuLista,TFLOAT >::operator =( destra); + VectorOp,TFLOAT >::operator =( destra); return *this; } -template unsigned int GreenKuboNComponentIonicFluid::numeroTimestepsOltreFineBlocco(unsigned int n_b) { +template unsigned int GreenKuboNComponentIonicFluid::nExtraTimesteps(unsigned int n_b) { return (traiettoria->n_timestep()/(n_b+1)+1 < lmax || lmax==0)? traiettoria->n_timestep()/(n_b+1)+1 : lmax; } template void GreenKuboNComponentIonicFluid::reset(unsigned int numeroTimestepsPerBlocco) { leff=(numeroTimestepsPerBlocco @@ -153,7 +153,7 @@ template unsigned int GreenKuboN return 3*N_corr; } -template void GreenKuboNComponentIonicFluid::calcola(unsigned int primo) { +template void GreenKuboNComponentIonicFluid::calculate(unsigned int primo) { if(!benchmarked) @@ -231,8 +231,8 @@ template void GreenKuboNComponen { //questa diventa la media della media nei pezzettini (devo aggiungere un contatore e usare la formula della media) //possibile perdita di precisione(?) - TFLOAT delta_JJ=JJ[j1]-lista[(itimestep)*narr+j1]; - lista[(itimestep)*narr+j1]+=delta_JJ/cont_JJ; + TFLOAT delta_JJ=JJ[j1]-vdata[(itimestep)*narr+j1]; + vdata[(itimestep)*narr+j1]+=delta_JJ/cont_JJ; } //N_corr (funzioni di correlazione), N_corr (integrali,integrali di einstein), 1 (kappa), 1 (kappa_einstein) @@ -261,7 +261,7 @@ template void GreenKuboNComponen for (unsigned int itimestep=start_mean;itimestep void GreenKuboNComponen if (subtract_mean) { //toglie la media a tutte le funzioni di correlazione prima di fare gli integrali for (unsigned int itimestep=0;itimestep void GreenKuboNComponen // I[i] = I[i-1]+ f[i-1]/2.0 + f[i]/2.0 for (unsigned int j=0;j void GreenKuboNComponen for (unsigned int j1=0;j1 void GreenKuboNComponen k= (coeff.block(0,0,1,1) - coeff.block(0,1,1,idx_j.size()-1)*coeff.block(1,1,idx_j.size()-1,idx_j.size()-1).inverse()*coeff.block(1,0,idx_j.size()-1,1))(0,0); else k=coeff(0,0); - lista[(itimestep)*narr+3*N_corr+0]=k; + vdata[(itimestep)*narr+3*N_corr+0]=k; //stessa cosa con la formula di einstein for (unsigned int j=0;j void GreenKuboNComponen for (unsigned int j1=0;j1 void GreenKuboNComponen k= (coeff.block(0,0,1,1) - coeff.block(0,1,1,idx_j.size()-1)*coeff.block(1,1,idx_j.size()-1,idx_j.size()-1).inverse()*coeff.block(1,0,idx_j.size()-1,1))(0,0); else k=coeff(0,0); - lista[(itimestep)*narr+3*N_corr+1]=k; + vdata[(itimestep)*narr+3*N_corr+1]=k; } @@ -351,7 +351,7 @@ template void GreenKuboNComponen //divide per itimestep tutti gli integrali einsteniani for (unsigned int itimestep=1;itimestep void GreenKuboNComponen #endif for (unsigned int itimestep=0;itimestep unsigned int GreenKuboN cronometro cron; cron.start(); - calcola(0); + calculate(0); cron.stop(); std::cerr << i<<" "< -#include - -HeatC::HeatC(Traiettoria *t, double sigma, unsigned int nthreads, unsigned int skip) : - sigma{sigma}, nthreads{nthreads},skip{skip},t{t},ntimesteps{0} -{ - if (nthreads==0) - nthreads=1; - if (skip==0) - skip=1; -} - -void HeatC::reset(const unsigned int numeroTimestepsPerBlocco) { - if (ntimesteps>= numeroTimestepsPerBlocco && t->get_ntypes()>=tipi_atomi){ - ntimesteps=numeroTimestepsPerBlocco; - return; - } - delete [] lista; - tipi_atomi=t->get_ntypes(); - ntimesteps=numeroTimestepsPerBlocco; - lunghezza_lista=ntimesteps*3*3*tipi_atomi/skip; - lista = new double[lunghezza_lista]; -} - -double calc0(const double & pre, const double & A,const double & B) { - return pre*(B*B-A*A); -} - -template -void calc(double * box,double * pos, const double * val, const double & sigma, double * res) { - const double sigma3=sigma*sigma*sigma*2; - for (unsigned int ival=0;ival<3;++ival){ - for (unsigned int icoord=0;icoord<3;++icoord) { - if (pos[icoord]+sigma/2.0>box[icoord*2+1]) { - res[3*ival+icoord]+=C*calc0(val[ival]/sigma3,pos[icoord]-sigma/2.0, box[icoord*2+1]); - res[3*ival+icoord]+=C*calc0(val[ival]/sigma3,box[icoord*2+0], pos[icoord]+sigma/2.0-box[icoord*2+1]+box[icoord*2+0]); - } else if (pos[icoord]-sigma/2.0 threads; - for (unsigned int ith=0;ithscatola(itimestep); - double *cellp = t->scatola(itimestep+1); - for (unsigned int iatom=0;iatomget_natoms();++iatom){ - double *rm = t->posizioni(itimestep,iatom); - double *rp = t->posizioni(itimestep+1,iatom); - double *vm = t->velocita(itimestep,iatom); - double *vp = t->velocita(itimestep+1,iatom); - calc<-1>(cellm,rm,vm,sigma,&lista[ilista*tipi_atomi*3*3+t->get_type(iatom)*3*3]); - calc<1>(cellp,rp,vp,sigma,&lista[ilista*tipi_atomi*3*3+t->get_type(iatom)*3*3]); - } - - ilista++; - } - })); - } - - for (auto & t : threads){ - t.join(); - } - threads.clear(); -} - -std::vector HeatC::get_shape() const{ - return {ntimesteps/skip,tipi_atomi,3,3}; -} -std::vector HeatC::get_stride() const{ - return {static_cast(sizeof(double)*tipi_atomi*3*3), - static_cast(sizeof(double)*3*3), - sizeof(double)*3, - sizeof(double)}; -} - -HeatC::~HeatC(){ - -} - -HeatC & HeatC::operator =(const HeatC & destra) { - OperazioniSuLista::operator =(destra); - return *this; -} diff --git a/lib/src/heatfluxts.cpp b/lib/src/heatfluxts.cpp index 1cd105d3c..96da592d7 100644 --- a/lib/src/heatfluxts.cpp +++ b/lib/src/heatfluxts.cpp @@ -17,7 +17,7 @@ //TODO: anche questo dovrebbe leggere il file a blocchi! -HeatFluxTs::HeatFluxTs(std::string filename, Traiettoria *t,unsigned int skip) +HeatFluxTs::HeatFluxTs(std::string filename, Trajectory *t,unsigned int skip) : traiettoria(t),heatflux(0),skip(skip) { diff --git a/lib/src/istogrammaatomiraggio.cpp b/lib/src/istogrammaatomiraggio.cpp index 7c4e2588d..ab37a3082 100644 --- a/lib/src/istogrammaatomiraggio.cpp +++ b/lib/src/istogrammaatomiraggio.cpp @@ -11,10 +11,10 @@ #include "istogrammaatomiraggio.h" -#include "traiettoria.h" +#include "trajectory.h" #include -IstogrammaAtomiRaggio::IstogrammaAtomiRaggio(Traiettoria *t, double r, unsigned int skip,unsigned int nthreads) : r2(r*r), skip(skip),traiettoria(t),nthreads(nthreads),hist(0) +IstogrammaAtomiRaggio::IstogrammaAtomiRaggio(Trajectory *t, double r, unsigned int skip,unsigned int nthreads) : r2(r*r), skip(skip),traiettoria(t),nthreads(nthreads),hist(0) { if (IstogrammaAtomiRaggio::skip<1) IstogrammaAtomiRaggio::skip=1; if (IstogrammaAtomiRaggio::nthreads<1) IstogrammaAtomiRaggio::nthreads=1; @@ -28,7 +28,7 @@ void IstogrammaAtomiRaggio::reset(const unsigned int numeroTimestepsPerBlocco) { hist= new std::map [ntypes]; } -void IstogrammaAtomiRaggio::calcola(unsigned int tstart) { +void IstogrammaAtomiRaggio::calculate(unsigned int tstart) { unsigned int natomith=natoms/nthreads; std::vector threads; diff --git a/lib/src/istogrammavelocita.cpp b/lib/src/istogrammavelocita.cpp index f10d9ec9e..008721d33 100644 --- a/lib/src/istogrammavelocita.cpp +++ b/lib/src/istogrammavelocita.cpp @@ -11,7 +11,6 @@ #include "istogrammavelocita.h" -#include "mediablocchi.h" void istg(double a,/// estremo inferiore dell'intervallo da istogrammare @@ -27,43 +26,43 @@ double k; } -IstogrammaVelocita::IstogrammaVelocita(Traiettoria *t, unsigned int nbins,double vminmax_) : OperazioniSuLista () +IstogrammaVelocita::IstogrammaVelocita(Trajectory *t, unsigned int nbins,double vminmax_) : VectorOp () { traiettoria=t; - lunghezza_lista=traiettoria->get_ntypes()*3*nbins; - lista = new double [lunghezza_lista]; + data_length=traiettoria->get_ntypes()*3*nbins; + vdata = new double [data_length]; bins=nbins; vminmax=vminmax_; } -unsigned int IstogrammaVelocita::numeroTimestepsOltreFineBlocco(unsigned int n_b) { +unsigned int IstogrammaVelocita::nExtraTimesteps(unsigned int n_b) { return 0; } IstogrammaVelocita & IstogrammaVelocita::operator = (const IstogrammaVelocita & destra) { - OperazioniSuLista::operator =(destra); + VectorOp::operator =(destra); return *this; } void IstogrammaVelocita::reset(const unsigned int numeroTimestepsPerBlocco) { - for (unsigned int i=0;iget_natoms();iatom++) { for (unsigned int icoord=0;icoord<3;icoord++) - istg(-vminmax,+vminmax,bins,&lista[(traiettoria->get_type(iatom)*3+icoord)*bins],traiettoria->velocita(itimestep,iatom)[icoord]); + istg(-vminmax,+vminmax,bins,&vdata[(traiettoria->get_type(iatom)*3+icoord)*bins],traiettoria->velocity(itimestep,iatom)[icoord]); } @@ -73,4 +72,4 @@ void IstogrammaVelocita::calcola(unsigned int primo) { -//template class MediaBlocchi; +//template class BlockAverage; diff --git a/lib/src/modivibrazionali.cpp b/lib/src/modivibrazionali.cpp index 047f7e5ce..359ae60cb 100644 --- a/lib/src/modivibrazionali.cpp +++ b/lib/src/modivibrazionali.cpp @@ -11,12 +11,12 @@ #include "modivibrazionali.h" -#include "traiettoria.h" #include #include #include #include #include +#include const double evjoules = 1.60217733e-19, // from eV to Joules @@ -37,7 +37,7 @@ void istg(double a,/// estremo inferiore dell'intervallo da istogrammare //unità immaginaria const std::complex I(0.0,1.0); -ModiVibrazionali::ModiVibrazionali(Traiettoria * tr, std::string ifcfile, std::string fononefile, unsigned int n_threads,unsigned int timestep_blocco) : OperazioniSuLista() +ModiVibrazionali::ModiVibrazionali(Trajectory * tr, std::string ifcfile, std::string fononefile, unsigned int n_threads,unsigned int timestep_blocco) : VectorOp() { if (n_threads==0) numero_threads=1; @@ -74,7 +74,7 @@ void ModiVibrazionali::read_force_file(std::string f) { abort(); } - // legge anche le masse degli atomi + // legge anche le mass degli atomi for (unsigned int i=0;iget_natoms()*3; - delete [] lista; - lista=new double[lunghezza_lista]; + data_length=traiettoria->get_natoms()*3; + delete [] vdata; + vdata=new double[data_length]; posizioni_equilibrio->reset(s); } void ModiVibrazionali::azzera() { - for (unsigned int i=0;icalcola(primo); +void ModiVibrazionali::calculate(unsigned int primo) { + posizioni_equilibrio->calculate(primo); /*calcola, per ogni atomo e per ogni modo normale, il prodotto scalare * fra lo spostamento dell'atomo dalla sua posizione d'equilibrio e @@ -236,7 +236,7 @@ void ModiVibrazionali::calcola(unsigned int primo) { } - // la radice delle masse è già inserica nel file ifc + // la radice delle mass è già inserica nel file ifc /* //moltiplica per il coefficiente 1/sqrt(m_i m_j) le comoponenti della matrice for (unsigned int iatom=0;iatomget_atoms_cell();iatom++) @@ -386,11 +386,11 @@ void ModiVibrazionali::calcola(unsigned int primo) { // se ho una traiettoria troppo lunga, che non sta nella ram in un colpo solo, devo caricare il file a pezzi (compito svolto da traiettoria) unsigned int ntimestep_load=timestepBlocco + numero_threads - timestepBlocco%numero_threads; if (timestepBlocco!=0) - traiettoria->imposta_dimensione_finestra_accesso(ntimestep_load); //quanti dati carico in memoria + traiettoria->set_data_access_block_size(ntimestep_load); //quanti dati carico in memoria for (unsigned int istep=0;istepimposta_inizio_accesso(istep+primo); // carica il prossimo blocco di dati + traiettoria->set_access_at(istep+primo); // carica il prossimo blocco di dati } @@ -436,7 +436,7 @@ void ModiVibrazionali::calcola(unsigned int primo) { q_punto += sqrt(traiettoria->get_mass( traiettoria->get_type(iatom) )) *exp(I*vettori_onda.col(imodo/(3*posizioni_equilibrio->get_atoms_cell())).dot( Eigen::Map(posizioni_equilibrio->get_fitted_pos(iatom)) - Eigen::Map(posizioni_equilibrio->get_atom_position_origin_cell(posizioni_equilibrio->get_atom_base_index(iatom))) ) )* (autovettori.block<3,1>(0+3*posizioni_equilibrio->get_atom_base_index(iatom),imodo).dot( - Eigen::Map(traiettoria->velocita(ith+primo,iatom))) ); + Eigen::Map(traiettoria->velocity(ith+primo,iatom))) ); // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ // velocity of iatom, current timestep } @@ -556,9 +556,9 @@ void ModiVibrazionali::calcola(unsigned int primo) { -/* questo probabilmente non serve, basta fare l'analisi a blocchi del risultato finale, senza includere anche le posizioni medie +/* questo probabilmente non serve, basta fare l'analisi a blocchi del risultato finale, senza includere anche le positions medie ModiVibrazionali & ModiVibrazionali::operator =(const ModiVibrazionali &d){ - OperazioniSuLista operator= (d); + VectorOp operator= (d); posizioni_equilibrio } diff --git a/lib/src/msd.cpp b/lib/src/msd.cpp index 60915d1ed..be988423f 100644 --- a/lib/src/msd.cpp +++ b/lib/src/msd.cpp @@ -17,7 +17,6 @@ #include #include "msd.h" #include "config.h" -#include "traiettoria.h" #include "floating_exceptions.h" #ifdef USE_MPI @@ -34,7 +33,7 @@ MSD::MSD(T *t, unsigned int skip, unsigned int tmax, unsigned int nthread f_cm=1; } template -unsigned int MSD::numeroTimestepsOltreFineBlocco(unsigned int n_b){ +unsigned int MSD::nExtraTimesteps(unsigned int n_b){ return (traiettoria->get_ntimesteps()/(n_b+1)+1 < lmax || lmax==0)? traiettoria->get_ntimesteps()/(n_b+1)+1 : lmax; } template @@ -42,11 +41,11 @@ void MSD::reset(const unsigned int numeroTimestepsPerBlocco) { leff=(numeroTimestepsPerBloccoget_ntypes(); - lunghezza_lista=leff*ntypes*f_cm; + data_length=leff*ntypes*f_cm; ntimesteps=numeroTimestepsPerBlocco; - delete [] lista; - lista=new double [lunghezza_lista]; + delete [] vdata; + vdata=new double [data_length]; } template @@ -67,7 +66,7 @@ void MSD::calc_single_th(size_t begin, size_t ultimo, size_t primo, unsig for (size_t t=begin;t::calc_single_th(size_t begin, size_t ultimo, size_t primo, unsig for (size_t iatom=0;iatomget_natoms();iatom++) { size_t itype=traiettoria->get_type(iatom); double delta=(pow( - traiettoria->template posizioni(primo+imedia,iatom)[0]-traiettoria->template posizioni(primo+imedia+t,iatom)[0] - -(traiettoria->template posizioni_cm(primo+imedia,itype)[0]-traiettoria->template posizioni_cm(primo+imedia+t,itype)[0]) + traiettoria->template positions(primo+imedia,iatom)[0]-traiettoria->template positions(primo+imedia+t,iatom)[0] + -(traiettoria->template positions_cm(primo+imedia,itype)[0]-traiettoria->template positions_cm(primo+imedia+t,itype)[0]) ,2)+ pow( - traiettoria->template posizioni(primo+imedia,iatom)[1]-traiettoria->template posizioni(primo+imedia+t,iatom)[1] - -(traiettoria->template posizioni_cm(primo+imedia,itype)[1]-traiettoria->template posizioni_cm(primo+imedia+t,itype)[1]) + traiettoria->template positions(primo+imedia,iatom)[1]-traiettoria->template positions(primo+imedia+t,iatom)[1] + -(traiettoria->template positions_cm(primo+imedia,itype)[1]-traiettoria->template positions_cm(primo+imedia+t,itype)[1]) ,2)+ pow( - traiettoria->template posizioni(primo+imedia,iatom)[2]-traiettoria->template posizioni(primo+imedia+t,iatom)[2] - -(traiettoria->template posizioni_cm(primo+imedia,itype)[2]-traiettoria->template posizioni_cm(primo+imedia+t,itype)[2]) + traiettoria->template positions(primo+imedia,iatom)[2]-traiettoria->template positions(primo+imedia+t,iatom)[2] + -(traiettoria->template positions_cm(primo+imedia,itype)[2]-traiettoria->template positions_cm(primo+imedia+t,itype)[2]) ,2)) - -lista[ntypes*t*f_cm+traiettoria->get_type(iatom)]; - lista[ntypes*t*f_cm+traiettoria->get_type(iatom)]+=delta/(++cont[traiettoria->get_type(iatom)]); + -vdata[ntypes*t*f_cm+traiettoria->get_type(iatom)]; + vdata[ntypes*t*f_cm+traiettoria->get_type(iatom)]+=delta/(++cont[traiettoria->get_type(iatom)]); } }else{ for (size_t iatom=0;iatomget_natoms();iatom++) { - double delta=(pow(traiettoria->template posizioni(primo+imedia,iatom)[0]-traiettoria->template posizioni(primo+imedia+t,iatom)[0],2)+ - pow(traiettoria->template posizioni(primo+imedia,iatom)[1]-traiettoria->template posizioni(primo+imedia+t,iatom)[1],2)+ - pow(traiettoria->template posizioni(primo+imedia,iatom)[2]-traiettoria->template posizioni(primo+imedia+t,iatom)[2],2)) - -lista[ntypes*t*f_cm+traiettoria->get_type(iatom)]; - lista[ntypes*t*f_cm+traiettoria->get_type(iatom)]+=delta/(++cont[traiettoria->get_type(iatom)]); - if constexpr (FPE) fpem.check_nan(lista[ntypes*t*f_cm+traiettoria->get_type(iatom)]); + double delta=(pow(traiettoria->template positions(primo+imedia,iatom)[0]-traiettoria->template positions(primo+imedia+t,iatom)[0],2)+ + pow(traiettoria->template positions(primo+imedia,iatom)[1]-traiettoria->template positions(primo+imedia+t,iatom)[1],2)+ + pow(traiettoria->template positions(primo+imedia,iatom)[2]-traiettoria->template positions(primo+imedia+t,iatom)[2],2)) + -vdata[ntypes*t*f_cm+traiettoria->get_type(iatom)]; + vdata[ntypes*t*f_cm+traiettoria->get_type(iatom)]+=delta/(++cont[traiettoria->get_type(iatom)]); + if constexpr (FPE) fpem.check_nan(vdata[ntypes*t*f_cm+traiettoria->get_type(iatom)]); } } if constexpr (FPE) { for (size_t itype=0;itypetemplate posizioni_cm(primo+imedia,itype)[0]-traiettoria->template posizioni_cm(primo+imedia+t,itype)[0],2)+ - pow(traiettoria->template posizioni_cm(primo+imedia,itype)[1]-traiettoria->template posizioni_cm(primo+imedia+t,itype)[1],2)+ - pow(traiettoria->template posizioni_cm(primo+imedia,itype)[2]-traiettoria->template posizioni_cm(primo+imedia+t,itype)[2],2)) - -lista[ntypes*t*f_cm+ntypes+itype]; - lista[ntypes*t*f_cm+ntypes+itype]+=delta/(++cont[ntypes+itype]); - if constexpr (FPE) fpem.check_nan(lista[ntypes*t*f_cm+ntypes+itype]); + double delta=(pow(traiettoria->template positions_cm(primo+imedia,itype)[0]-traiettoria->template positions_cm(primo+imedia+t,itype)[0],2)+ + pow(traiettoria->template positions_cm(primo+imedia,itype)[1]-traiettoria->template positions_cm(primo+imedia+t,itype)[1],2)+ + pow(traiettoria->template positions_cm(primo+imedia,itype)[2]-traiettoria->template positions_cm(primo+imedia+t,itype)[2],2)) + -vdata[ntypes*t*f_cm+ntypes+itype]; + vdata[ntypes*t*f_cm+ntypes+itype]+=delta/(++cont[ntypes+itype]); + if constexpr (FPE) fpem.check_nan(vdata[ntypes*t*f_cm+ntypes+itype]); } } } @@ -138,7 +137,7 @@ void MSD::calc_end() { for (size_t ts=0;ts::calc_end() { template MSD & MSD::operator=(const MSD &destra) { - OperazioniSuLista >::operator =( destra); + VectorOp >::operator =( destra); return *this; } #ifdef BUILD_MMAP -template class MSD; -template class MSD; +#include "trajectory.h" +template class MSD; +template class MSD; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class MSD; -template class MSD; +#include "trajectory_numpy.h" +template class MSD; +template class MSD; #endif diff --git a/lib/src/neighbour.cpp b/lib/src/neighbour.cpp index 675b06453..200d6c1bc 100644 --- a/lib/src/neighbour.cpp +++ b/lib/src/neighbour.cpp @@ -99,10 +99,10 @@ size_t Neighbours::get_sann_n(const size_t iatom, const size_t jtype) c } #ifdef BUILD_MMAP -#include "traiettoria.h" -template class Neighbours; +#include "trajectory.h" +template class Neighbours; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class Neighbours; +#include "trajectory_numpy.h" +template class Neighbours; #endif diff --git a/lib/src/operazionisulista.cpp b/lib/src/operazionisulista.cpp index ddd6482fb..cd4406fbe 100644 --- a/lib/src/operazionisulista.cpp +++ b/lib/src/operazionisulista.cpp @@ -22,51 +22,45 @@ #include "greenkuboNcomponentionicfluid.h" #include "gofrt.h" #include "msd.h" -#include "correlatorespaziale.h" -#include "heatc.h" -#include "calcolamultithread.h" -#include "centerdiff.h" -#include "centerofmassdiff.h" - -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator + (const T & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator + (const T & destra) const { T res = static_cast(*this); res+=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator - (const T & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator - (const T & destra) const { T res = static_cast(*this); res-=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator * (const T & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator * (const T & destra) const { T res = static_cast(*this); res*=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator / (const T & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator / (const T & destra) const { T res = static_cast(*this); res/=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator + (const TFLOAT & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator + (const TFLOAT & destra) const { T res = static_cast(*this); res+=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator - (const TFLOAT & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator - (const TFLOAT & destra) const { T res = static_cast(*this); res-=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator * (const TFLOAT & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator * (const TFLOAT & destra) const { T res = static_cast(*this); res*=destra; return res; } -template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::operator / (const TFLOAT & destra) const { +template < class T, class TFLOAT > const T VectorOp < T, TFLOAT >::operator / (const TFLOAT & destra) const { T res = static_cast(*this); res/=destra; return res; @@ -75,20 +69,16 @@ template < class T, class TFLOAT > const T OperazioniSuLista < T, TFLOAT >::ope -template class OperazioniSuLista; -template class OperazioniSuLista; -template class OperazioniSuLista; -template class OperazioniSuLista; -template class OperazioniSuLista; -template class OperazioniSuLista >; -template class OperazioniSuLista,long double >; -template class OperazioniSuLista,long double >; -template class OperazioniSuLista; -template class OperazioniSuLista,double>; -template class OperazioniSuLista,long double>; -template class OperazioniSuLista; -template class OperazioniSuLista; -template class OperazioniSuLista; -template class OperazioniSuLista; -//template class OperazioniSuLista; +template class VectorOp; +template class VectorOp; +template class VectorOp; +template class VectorOp; +template class VectorOp; +template class VectorOp >; +template class VectorOp,long double >; +template class VectorOp,long double >; +template class VectorOp; +template class VectorOp,double>; +template class VectorOp,long double>; +template class VectorOp; */ diff --git a/lib/src/posizioniequilibrio.cpp b/lib/src/posizioniequilibrio.cpp index 9b76f3142..6bb2ed3fd 100644 --- a/lib/src/posizioniequilibrio.cpp +++ b/lib/src/posizioniequilibrio.cpp @@ -11,22 +11,21 @@ #include "posizioniequilibrio.h" -#include "traiettoria.h" #include #include #include #include #include "doubleround.h" -PosizioniEquilibrio::PosizioniEquilibrio(Traiettoria * tr,unsigned int timesteps_sottoblocco) +PosizioniEquilibrio::PosizioniEquilibrio(Trajectory * tr,unsigned int timesteps_sottoblocco) { traiettoria=tr; - lunghezza_lista=traiettoria->get_natoms()*3; - lista = new double [lunghezza_lista]; - posizioni_fittate_reticolo=new double[lunghezza_lista]; - for (unsigned int i=0;iget_natoms()*3; + vdata = new double [data_length]; + posizioni_fittate_reticolo=new double[data_length]; + for (unsigned int i=0;i > [traiettoria->get_natoms()]; calcolato=false; @@ -103,7 +102,7 @@ PosizioniEquilibrio::~PosizioniEquilibrio(){ delete [] traslation; } -unsigned int PosizioniEquilibrio::numeroTimestepsOltreFineBlocco(unsigned int n_b){ +unsigned int PosizioniEquilibrio::nExtraTimesteps(unsigned int n_b){ return 0; } @@ -111,23 +110,23 @@ void PosizioniEquilibrio::reset(const unsigned int numeroTimestepsPerBlocco){ lunghezza_media=numeroTimestepsPerBlocco; } -void PosizioniEquilibrio::calcola(unsigned int primo) { +void PosizioniEquilibrio::calculate(unsigned int primo) { azzera(); if (timestepSottoblocco>0) - traiettoria->imposta_dimensione_finestra_accesso(timestepSottoblocco); + traiettoria->set_data_access_block_size(timestepSottoblocco); for (unsigned int i=0;i0 && i%timestepSottoblocco==0){ - traiettoria->imposta_inizio_accesso(primo+i); + traiettoria->set_access_at(primo+i); } for (unsigned int iatom=0;iatomget_natoms();iatom++){ for (unsigned int icoord=0;icoord<3;icoord++){ - double delta = traiettoria->posizioni(primo+i,iatom)[icoord] - - lista[iatom*3+icoord]; - lista[iatom*3+icoord] += delta/(i+1); + double delta = traiettoria->positions(primo+i,iatom)[icoord] + - vdata[iatom*3+icoord]; + vdata[iatom*3+icoord] += delta/(i+1); } } } @@ -141,7 +140,7 @@ double * PosizioniEquilibrio::get_atom_position(unsigned int iatom){ return 0; } - return & lista[iatom*3]; + return & vdata[iatom*3]; } void PosizioniEquilibrio::reticolo_inizializza(double *base_r, double *base,unsigned int *base_type, const unsigned int &nbase) { @@ -289,13 +288,13 @@ void PosizioniEquilibrio::coord_reticolo(double *xyz, double *uvw_min, double *u -double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double *max, /// limiti della scatola +double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double *max, /// limiti della box double *spostamento ) { //obiettivo: voglio assegnare ad ogni atomo il suo indice di cella, - //e calcolare esattamente le posizioni degli atomi secondo il reticolo a temperatura zero + //e calcolare esattamente le positions degli atomi secondo il reticolo a temperatura zero - //assumo che il reticolo sia fisso. devo solo trovare come traslarlo per centrare nel migliore dei modi gli atomi, che non si trovano esattamente sulle posizioni di temperatura zero (la temperatura della simulazione è finita) + //assumo che il reticolo sia fisso. devo solo trovare come traslarlo per centrare nel migliore dei modi gli atomi, che non si trovano esattamente sulle positions di temperatura zero (la temperatura della simulazione è finita) //sceglie un atomo di uno spigolo e lo considera come origine. devo controllare di poter costruire una cella con questo atomo //gli atomi dello spigolo sono quelli che hanno meno primi vicini. @@ -311,7 +310,7 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * if (iatom!=jatom) { double d2=0.0; for (unsigned int i=0;i<3;i++) - d2+=(lista[iatom*3+i]-lista[jatom*3+i])*(lista[iatom*3+i]-lista[jatom*3+i]); + d2+=(vdata[iatom*3+i]-vdata[jatom*3+i])*(vdata[iatom*3+i]-vdata[jatom*3+i]); distanze.insert(std::pair(d2,jatom)); } } @@ -343,8 +342,8 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * for (unsigned int iatom=0;iatomget_natoms();iatom++) { double d2=0.0; for (unsigned int i=0;i<3;i++ ) - d2+=(lista[iatom*3+i]-(lista[*iterator *3+i]+base_reticolo(i,ibase)) )* - (lista[iatom*3+i]-(lista[*iterator *3+i]+base_reticolo(i,ibase)) ); + d2+=(vdata[iatom*3+i]-(vdata[*iterator *3+i]+base_reticolo(i,ibase)) )* + (vdata[iatom*3+i]-(vdata[*iterator *3+i]+base_reticolo(i,ibase)) ); distanze.insert(std::pair(d2,iatom)); } //prendi il più vicino. se la distanza è maggiore di una cerca soglia, scarta questo atomo dell'angolo @@ -371,7 +370,7 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * abort(); } - double origin[3]={lista[3* *iterator+0],lista[3* *iterator+1],lista[3* *iterator+2]}; + double origin[3]={vdata[3* *iterator+0],vdata[3* *iterator+1],vdata[3* *iterator+2]}; //stima molto rozza dei limiti degli indici che coprono completamente la cella di simulazione double coord[]={max[0]-origin[0],max[1]-origin[1],max[2]-origin[2]}, @@ -397,7 +396,7 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * coord_reticolo(coord,uvw_min,uvw_max); // m,M,M /*adesso ho i limiti nelle nuove coordinate, genero i punti del reticolo - * nella scatola e poi calcolo la distanza da quelli della traiettoria + * nella box e poi calcolo la distanza da quelli della traiettoria */ spostamento[0]=0.0; @@ -409,7 +408,7 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * for (int u=std::lrint( uvw_min[0]);u<=std::lrint(uvw_max[0]);u++){ for (int v=std::lrint( uvw_min[1]);v<=std::lrint(uvw_max[1]);v++){ for (int w=std::lrint( uvw_min[2]);w<=std::lrint(uvw_max[2]);w++){ - //controlla che l'atomo sia dentro la scatola, se si lo aggiunge + //controlla che l'atomo sia dentro la box, se si lo aggiunge //per ogni atomo della base... Eigen::Vector3d uvw; uvw << double(u),double(v),double(w); @@ -430,7 +429,7 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * } if (reticolo_xyz.size()< traiettoria->get_natoms()) { - std::cerr << "La scatola non è abbastanza grande per contenere il reticolo intero ("<get_natoms()<<" atomi della traiettoria)!\n"; + std::cerr << "La box non è abbastanza grande per contenere il reticolo intero ("<get_natoms()<<" atomi della traiettoria)!\n"; } //guardo la distanza dei punti da quelli della media (da calcolare in precedenza) @@ -439,12 +438,12 @@ double PosizioniEquilibrio::d2_reticolo_spostamento_medio(double * min, double * std::array uvwi=reticolo_xyz.begin()->first; // cerco l'atomo più vicino dello stesso tipo, per restituire comunque un numero sensato auto iteratore = reticolo_xyz.begin(); - Eigen::Vector3d d = iteratore->second.first - Eigen::Map(&lista[iatom*3]); + Eigen::Vector3d d = iteratore->second.first - Eigen::Map(&vdata[iatom*3]); Eigen::Vector3d dmin; d2m= d.dot(d)*4; for (iteratore++;iteratore != reticolo_xyz.end();iteratore++) { if (base_tipi[iteratore->first[3]]==traiettoria->get_type(iatom) && iteratore->second.second == traiettoria->get_natoms()){ - d = iteratore->second.first - Eigen::Map(&lista[iatom*3]); + d = iteratore->second.first - Eigen::Map(&vdata[iatom*3]); d2t= d.dot(d); if (d2tget_natoms()) { - std::cerr << "Errore: questo punto del reticolo ha già un atomo assegnato (provare ad aumentare le dimensioni della scatola dove si genera il reticolo?)\n"; + std::cerr << "Errore: questo punto del reticolo ha già un atomo assegnato (provare ad aumentare le dimensioni della box dove si genera il reticolo?)\n"; abort(); } reticolo_xyz[uvwi].second=iatom; @@ -497,7 +496,7 @@ double * PosizioniEquilibrio::get_fitted_pos(unsigned int iatom) { void PosizioniEquilibrio::get_displacement(unsigned int iatom, unsigned int tstep,double * displ) { for (unsigned int i=0;i<3;i++) { - displ[i]=traiettoria->posizioni(tstep,iatom)[i]-get_fitted_pos(iatom)[i]; + displ[i]=traiettoria->positions(tstep,iatom)[i]-get_fitted_pos(iatom)[i]; } } @@ -557,7 +556,7 @@ unsigned int PosizioniEquilibrio::get_number_cells(){ double PosizioniEquilibrio::get_simulation_size(unsigned int icoord) { if (icoord <3) { - double *b=traiettoria->scatola_last(); + double *b=traiettoria->box_last(); return b[icoord*2+1]-b[icoord*2]; } else { std::cerr << "Errore: richiesta una dimensione che non esiste in get_simulation_size()"; @@ -579,13 +578,13 @@ void PosizioniEquilibrio::fit_nacl(){ * Lo spostamento di prima più questo nuovo sarà lo spostamento effettivo del reticolo rispetto all'origine. */ - //trova le dimensioni della scatola, leggendole dalla traiettoria - double max[3]={traiettoria->scatola_last()[1],traiettoria->scatola_last()[3],traiettoria->scatola_last()[5]}, - min[3]={traiettoria->scatola_last()[0],traiettoria->scatola_last()[2],traiettoria->scatola_last()[4]}; + //trova le dimensioni della box, leggendole dalla traiettoria + double max[3]={traiettoria->box_last()[1],traiettoria->box_last()[3],traiettoria->box_last()[5]}, + min[3]={traiettoria->box_last()[0],traiettoria->box_last()[2],traiettoria->box_last()[4]}; #ifdef DEBUG - std::cerr << "Limiti della scatola (minxyz),(maxxyz) ("<< min[0]<< ", " << min[1] << ", " << min[2] << "), (" << max[0] << ", " << max[1] << ", " << max[2]<<")\n"; + std::cerr << "Limiti della box (minxyz),(maxxyz) ("<< min[0]<< ", " << min[1] << ", " << min[2] << "), (" << max[0] << ", " << max[1] << ", " << max[2]<<")\n"; #endif /* ______ @@ -596,7 +595,7 @@ void PosizioniEquilibrio::fit_nacl(){ * */ - //assumo scatola quadrata + //assumo box quadrata double media_primi_vicini= cbrt( (max[0]-min[0])*(max[1]-min[1])*(max[2]-min[2])/(traiettoria->get_natoms()/8))/2.0; for (unsigned int i=0;i<3;i++) { diff --git a/lib/src/readlog.cpp b/lib/src/readlog.cpp index b894d7401..3fb09bfb9 100644 --- a/lib/src/readlog.cpp +++ b/lib/src/readlog.cpp @@ -23,7 +23,7 @@ #include "chargefluxts.h" -template ReadLog::ReadLog(std::string filename, Traiettoria *t, unsigned int skip, unsigned int nthreads, unsigned int nbatch, std::vector req_headers): +template ReadLog::ReadLog(std::string filename, Trajectory *t, unsigned int skip, unsigned int nthreads, unsigned int nbatch, std::vector req_headers): traiettoria(t),skip(skip),nthreads(nthreads),nbatch(nbatch) { cronometro cron; @@ -84,7 +84,7 @@ template ReadLog::ReadLog(std::string filename, Traiettor } #endif - //determina il numero di colonne da allocare, e comprende le eventuali cariche + //determina il numero di colonne da allocare, e comprende le eventuali charge unsigned int n_columns_from_binary=0; data_size_from_binary=0; if (req_headers.size()>0){ @@ -94,7 +94,7 @@ template ReadLog::ReadLog(std::string filename, Traiettor } auto Qs=qs(*it); if (Qs.first !=""){ - n_columns_from_binary++; //aggiungi i valori delle cariche letti + n_columns_from_binary++; //aggiungi i valori delle charge letti q_current_type.push_back(Qs); } } @@ -220,7 +220,7 @@ template unsigned int ReadLog::get_calc_j_index(std::stri } -//questo analizza la stringa speciale "#traj:JZ N q1 ... qN" e ritorna le cariche +//questo analizza la stringa speciale "#traj:JZ N q1 ... qN" e ritorna le charge template std::pair > ReadLog::qs(std::string header) { if (header.size()==0 || header[0] != '#') return std::pair >(); @@ -257,7 +257,7 @@ template int ReadLog::need_binary(std::vector void ReadLog::calc_currents(Traiettoria * t,unsigned int n_b){ +template void ReadLog::calc_currents(Trajectory * t,unsigned int n_b){ traiettoria=t; //calcola e legge la corrente partendo dal file binario (nella classe traiettoria sono già presenti le velocità dei centri di massa) unsigned int timesteps_tot=data.size()/data_size; @@ -269,18 +269,18 @@ template void ReadLog::calc_currents(Traiettoria * t,unsi throw std::runtime_error("Cannot calculate current: number of coefficients must be equal to the number of atomic types in the trajectory!\n"); } } - traiettoria->imposta_dimensione_finestra_accesso(n_data_b); + traiettoria->set_data_access_block_size(n_data_b); for (unsigned int ib=0;ibimposta_dimensione_finestra_accesso(ultimo-ib*n_data_b); + traiettoria->set_data_access_block_size(ultimo-ib*n_data_b); } - traiettoria->imposta_inizio_accesso(n_data_b*ib); + traiettoria->set_access_at(n_data_b*ib); for (unsigned int ts=ib*n_data_b;tsvelocita_cm(ts,0); + double * v_cm=t->velocity_cm(ts,0); for (unsigned int icoord=0;icoord<3;icoord++) data[ts*data_size+(data_size-data_size_from_binary)+i*3+icoord]=0.0; for (unsigned int icm=0;icm +#include template -SpettroVibrazionale::SpettroVibrazionale(T * t, bool dump_):OperazioniSuLista() +SpettroVibrazionale::SpettroVibrazionale(T * t, bool dump_):VectorOp() { traiettoria=t; dump=dump_; @@ -48,28 +48,24 @@ void SpettroVibrazionale::deallocate_plan(){ } template -unsigned int SpettroVibrazionale::numeroTimestepsOltreFineBlocco(unsigned int n_b){ +unsigned int SpettroVibrazionale::nExtraTimesteps(unsigned int n_b){ return 0; } template void SpettroVibrazionale::reset(const unsigned int numeroTimestepsPerBlocco) { //inizializzo la memoria per i moduli quadri e basta, se necessario! -// size è la lunghezza della trasformata. La lista che contiene i moduli quadri sarà size/2+1 (trasformata reale). +// size è la lunghezza della trasformata. La vdata che contiene i moduli quadri sarà size/2+1 (trasformata reale). if (numeroTimestepsPerBlocco!=size) { size=numeroTimestepsPerBlocco; - delete [] lista; + delete [] vdata; tipi_atomi=traiettoria->get_ntypes(); if (tipi_atomi<=0) { + throw std::runtime_error("Number of types in the trajectory must be > 0"); tipi_atomi=1; - std::cerr << "Attenzione: non ho letto il numero di tipi diversi di atomi (non hai caricato la traiettoria prima di iniziare l'analisi?\n"; } - lunghezza_lista=(size/2+1)*3*tipi_atomi; // uno per ogni direzione dello spazio, per testare l'isotropia, e per tipo di atomo - lista = new double[lunghezza_lista]; -#ifdef DEBUG - std::cerr << "called SpettroVibrazionale::reset " __FILE__ ":"<<__LINE__<<"\n"; - std::cerr << "new double [] "<::reset(const unsigned int numeroTimestepsPerBlocco) //prima di chiamare questa la traiettoria deve essere impostata correttamente sulla finestra giusta (la funzione non sa qual'è lo timestep corrente. template -void SpettroVibrazionale::calcola(unsigned int primo ///ignorato: prendo l'inizio di quello che c'è in memoria (attenzione a caricare bene i dati!) +void SpettroVibrazionale::calculate(unsigned int primo ///ignorato: prendo l'inizio di quello che c'è in memoria (attenzione a caricare bene i dati!) ) { //alloca se necessario con fftw_malloc (che alloca la memoria in modo che sia allineata correttaemente per poter sfruttare le istruzioni SIMD del processore // e prepara il piano della trasformata. Trasformata_size è la dimensione della trasformata. Deve essere uguale a size, la dimensione dell'array dei moduli quadri. @@ -92,7 +88,7 @@ void SpettroVibrazionale::calcola(unsigned int primo ///ignorato: prendo l'i (int*)&size, // lunghezza di ciascuna trasformata 3*traiettoria->get_natoms(), // numero di trasformate /* ****** input ****** */ - dummy,//traiettoria->velocita_inizio(), // puntatore ai dati + dummy,//traiettoria->velocity_data(), // puntatore ai dati NULL, // i dati sono tutti compatti, non fanno parte di array più grandi 3*traiettoria->get_natoms(), // la distanza fra due dati successivi 1, // la distanza fra due serie di dati adiacenti @@ -109,13 +105,13 @@ void SpettroVibrazionale::calcola(unsigned int primo ///ignorato: prendo l'i //devo fare la trasformata della velocità per ogni atomo } - fftw_execute_dft_r2c(fftw3,traiettoria->velocita_inizio(),trasformata); + fftw_execute_dft_r2c(fftw3,traiettoria->velocity_data(),trasformata); // adesso bisogna fare la media del modulo quadro, una per ogni tipo di atomo double * media = new double[tipi_atomi]; unsigned int * cont = new unsigned int [tipi_atomi]; - for (unsigned int iomega=0;iomega::calcola(unsigned int primo ///ignorato: prendo l'i double modulo2=x*x+y*y; int itipo=traiettoria->get_type(iatom); if (itipo < 0 ) { - std::cerr << "Errore: tipo dell'atomo fuori dal range! (memoria corrotta?)\n"; - abort(); + throw std::runtime_error("Type index of the atom out of range (!?)"); itipo=0; } else if (itipo >= tipi_atomi) { - std::cerr << "Errore: tipo dell'atomo fuori dal range! (numerazione dei tipi non consecutiva?)\n"; - abort(); + throw std::runtime_error("Type index of the atom out of range (non consecutive indexes?!)"); itipo=tipi_atomi-1; } double delta = modulo2-media[itipo]; @@ -151,7 +145,7 @@ void SpettroVibrazionale::calcola(unsigned int primo ///ignorato: prendo l'i } for (unsigned int itipo=0;itipo < tipi_atomi;itipo++){ - lista[itipo*lunghezza_lista/tipi_atomi+ iomega*3+idim]=media[itipo]; + vdata[itipo*data_length/tipi_atomi+ iomega*3+idim]=media[itipo]; } } } @@ -160,7 +154,7 @@ void SpettroVibrazionale::calcola(unsigned int primo ///ignorato: prendo l'i if (dump) { std::ofstream out("analisi_vibr.debug",std::fstream::app); - for (unsigned int iomega=0;iomega::calcola(unsigned int primo ///ignorato: prendo l'i template double SpettroVibrazionale::spettro(unsigned int frequenza, unsigned int dim,unsigned int tipo_atomo) { - if (frequenza::spettro(unsigned int frequenza, unsigned int dim, template SpettroVibrazionale & SpettroVibrazionale::operator = (const SpettroVibrazionale & destra) { -#ifdef DEBUG - std::cerr << "called SpettroVibrazionale::operator =" __FILE__ ":"<<__LINE__<<"\n"; -#endif - OperazioniSuLista>::operator =(destra); + VectorOp>::operator =(destra); - //TODO: cos'è sta roba????? fftw_free(trasformata); // deallocate_plan(); trasformata=0; @@ -196,10 +186,11 @@ SpettroVibrazionale & SpettroVibrazionale::operator = (const SpettroVibraz } #ifdef BUILD_MMAP -template class SpettroVibrazionale; +#include "trajectory.h" +template class SpettroVibrazionale; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class SpettroVibrazionale; +#include "trajectory_numpy.h" +template class SpettroVibrazionale; #endif diff --git a/lib/src/sphericalbase.cpp b/lib/src/sphericalbase.cpp index 54154bbe2..1fa65eab6 100644 --- a/lib/src/sphericalbase.cpp +++ b/lib/src/sphericalbase.cpp @@ -99,14 +99,14 @@ void SphericalBase::calc(int timestep, } #include "config.h" #ifdef BUILD_MMAP -#include "traiettoria.h" -template class SphericalBase<6,double,Traiettoria>; -template class SphericalBase<8,double,Traiettoria>; -template class SphericalBase<10,double,Traiettoria>; +#include "trajectory.h" +template class SphericalBase<6,double,Trajectory>; +template class SphericalBase<8,double,Trajectory>; +template class SphericalBase<10,double,Trajectory>; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class SphericalBase<6,double,Traiettoria_numpy>; -template class SphericalBase<8,double,Traiettoria_numpy>; -template class SphericalBase<10,double,Traiettoria_numpy>; +#include "trajectory_numpy.h" +template class SphericalBase<6,double,Trajectory_numpy>; +template class SphericalBase<8,double,Trajectory_numpy>; +template class SphericalBase<10,double,Trajectory_numpy>; #endif diff --git a/lib/src/sphericalcorrelations.cpp b/lib/src/sphericalcorrelations.cpp index 5be716ef0..1ed987b3f 100644 --- a/lib/src/sphericalcorrelations.cpp +++ b/lib/src/sphericalcorrelations.cpp @@ -30,7 +30,7 @@ SphericalCorrelations::~SphericalCorrelations() { template -unsigned int SphericalCorrelations::numeroTimestepsOltreFineBlocco(unsigned int n_b) { +unsigned int SphericalCorrelations::nExtraTimesteps(unsigned int n_b) { return (t.get_ntimesteps()/(n_b+1)+1 < tmax || tmax==0)? t.get_ntimesteps()/(n_b+1)+1 : tmax; } @@ -49,18 +49,18 @@ void SphericalCorrelations::reset(const unsigned int numeroTimesteps //quanti timestep è lunga la funzione di correlazione leff =(numeroTimestepsPerBlocco -void SphericalCorrelations::calcola(unsigned int primo) { +void SphericalCorrelations::calculate(unsigned int primo) { if (leff+ntimesteps+primo > t.get_ntimesteps()){ throw std::runtime_error("trajectory is too short for this kind of calculation. Select a different starting timestep or lower the size of the average or the lenght of the correlation function"); @@ -98,7 +98,7 @@ void SphericalCorrelations::calcola(unsigned int primo) { TwoLoopSplit task_distributer(nthreads,ntimesteps,skip,skip*10,leff,1,block_t); //allocate space for per thread averages - TFLOAT * lista_th = new TFLOAT[lunghezza_lista*(nthreads-1)]; + TFLOAT * lista_th = new TFLOAT[data_length*(nthreads-1)]; unsigned int * lista_th_counters = new unsigned int [leff*nthreads]; size_t * hit = new size_t[nthreads]; size_t * miss = new size_t[nthreads]; @@ -110,10 +110,10 @@ void SphericalCorrelations::calcola(unsigned int primo) { if (neighList.size()>0) { nns = new Neighbours_T{&t,neighList}; } - TFLOAT * lista_th_= ith >0 ? lista_th+lunghezza_lista*(ith-1) : lista; + TFLOAT * lista_th_= ith >0 ? lista_th+data_length*(ith-1) : vdata; unsigned int * lista_th_counters_ = lista_th_counters+leff*ith; - for (unsigned int i=0;i::calcola(unsigned int primo) { //the data of the first thread is in place for (unsigned int dt=0;dt::corr_sh_calc(const TFLOAT * sh1, cons } #ifdef BUILD_MMAP -template class SphericalCorrelations<10,double,Traiettoria>; -template class SphericalCorrelations<6,double,Traiettoria>; +#include "trajectory.h" +template class SphericalCorrelations<10,double,Trajectory>; +template class SphericalCorrelations<6,double,Trajectory>; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class SphericalCorrelations<10,double,Traiettoria_numpy>; -template class SphericalCorrelations<6,double,Traiettoria_numpy>; +#include "trajectory_numpy.h" +template class SphericalCorrelations<10,double,Trajectory_numpy>; +template class SphericalCorrelations<6,double,Trajectory_numpy>; #endif diff --git a/lib/src/steinhardt.cpp b/lib/src/steinhardt.cpp index 423fe5552..d110ed67b 100644 --- a/lib/src/steinhardt.cpp +++ b/lib/src/steinhardt.cpp @@ -1,7 +1,7 @@ #include "steinhardt.h" #include "specialfunctions.h" - #include "config.h" +#include template Steinhardt::Steinhardt(T *t, @@ -15,7 +15,7 @@ Steinhardt::Steinhardt(T *t, const NeighListSpec nls, const bool averaged_order) : SPB::SphericalBase{t,nbin,rminmax}, - CMT::CalcolaMultiThread{nthreads,skip}, + CMT::CalculateMultiThread{nthreads,skip}, steinhardt_histogram_size{0}, steinhardt_l_histogram{steinhardt_l_histogram}, nbin_steinhardt{nbin_steinhardt}, @@ -39,7 +39,7 @@ Steinhardt::Steinhardt(T *t, for (unsigned i=0;i::reset(const unsigned int numeroTimestepsPerBlocco) c_descr=descr.str(); - if (new_lungh_lista != lunghezza_lista) { - delete [] lista; - lunghezza_lista=new_lungh_lista; - lista = new TFLOAT [lunghezza_lista]; + if (new_lungh_lista != data_length) { + delete [] vdata; + data_length=new_lungh_lista; + vdata = new TFLOAT [data_length]; } } @@ -87,7 +87,7 @@ void Steinhardt::calc_init(int primo) { //allocate the space for the threads work if (do_histogram){ - threadResults = new TFLOAT [lunghezza_lista*(CMT::nthreads-1)]; + threadResults = new TFLOAT [data_length*(CMT::nthreads-1)]; } else { threadResults = nullptr; } @@ -108,17 +108,17 @@ void Steinhardt::calc_single_th(int istart,//timestep index, begin TFLOAT *result = new TFLOAT [(l+1)*(l+1)*natoms*nbin*ntypes]; TFLOAT *result_averaged = nullptr; if (averaged_order) result_averaged = new TFLOAT [(l+1)*(l+1)*natoms*nbin*ntypes]; - TFLOAT * threadResult = threadResults + lunghezza_lista*ith; + TFLOAT * threadResult = threadResults + data_length*ith; int * counter = new int[natoms*ntypes*nbin]; if (do_histogram){ if (ith==CMT::nthreads-1) { // last thread writes directly on the final result memory - threadResult = lista; + threadResult = vdata; } - for (size_t i=0;i::join_data() { if (do_histogram){ //sum all the threads result for (int ith=0;ith::join_data() { } #ifdef BUILD_MMAP -template class Steinhardt<6,double,Traiettoria>; -template class Steinhardt<8,double,Traiettoria>; -template class Steinhardt<10,double,Traiettoria>; +#include "trajectory.h" +template class Steinhardt<6,double,Trajectory>; +template class Steinhardt<8,double,Trajectory>; +template class Steinhardt<10,double,Trajectory>; #endif #ifdef PYTHON_SUPPORT -#include "traiettoria_numpy.h" -template class Steinhardt<6,double,Traiettoria_numpy>; -template class Steinhardt<8,double,Traiettoria_numpy>; -template class Steinhardt<10,double,Traiettoria_numpy>; +#include "trajectory_numpy.h" +template class Steinhardt<6,double,Trajectory_numpy>; +template class Steinhardt<8,double,Trajectory_numpy>; +template class Steinhardt<10,double,Trajectory_numpy>; #endif diff --git a/lib/src/testtraiettoria.cpp b/lib/src/testtraiettoria.cpp index 06515b9ce..acc481cf6 100644 --- a/lib/src/testtraiettoria.cpp +++ b/lib/src/testtraiettoria.cpp @@ -13,7 +13,7 @@ #include "testtraiettoria.h" #include #include -TestTraiettoria::TestTraiettoria(std::string filename) : Traiettoria(filename) +TestTraiettoria::TestTraiettoria(std::string filename) : Trajectory(filename) { std::ofstream out("analisi_test.debug"); @@ -31,14 +31,14 @@ TestTraiettoria::TestTraiettoria(std::string filename) : Traiettoria(filename) unsigned int ts_b=n_timesteps/(n_b+1); - imposta_dimensione_finestra_accesso(ts_b*2); + set_data_access_block_size(ts_b*2); for (unsigned int ib=0;ib #include #include @@ -41,13 +41,13 @@ -Traiettoria::Traiettoria(std::string filename) +Trajectory::Trajectory(std::string filename) { struct stat sb; fsize=0; file=0; - timestep_corrente=0; + current_timestep=0; timestep_indicizzato=0; loaded_timesteps=0; tstep_size=0; @@ -55,12 +55,12 @@ Traiettoria::Traiettoria(std::string filename) n_timesteps=0; timesteps=0; timesteps_lammps=0; - buffer_tipi=0; - buffer_tipi_id=0; - buffer_posizioni=0; - buffer_velocita=0; - buffer_scatola=0; - buffer_posizioni_size=0; + buffer_type=0; + buffer_type_id=0; + buffer_positions=0; + buffer_velocity=0; + buffer_boxes=0; + buffer_positions_size=0; buffer_cm_size=0; ntypes=0; ok=false; @@ -68,12 +68,12 @@ Traiettoria::Traiettoria(std::string filename) indexed_all=false; min_type=0; max_type=0; - masse=0; - cariche=0; + mass=0; + charge=0; wrap_pbc=false; calculate_center_of_mass=false; - buffer_posizioni_cm=0; - buffer_velocita_cm=0; + buffer_positions_cm=0; + buffer_velocity_cm=0; fd=open(filename.c_str(), O_RDONLY); if (fd==-1) { @@ -103,11 +103,11 @@ Traiettoria::Traiettoria(std::string filename) size_t dimensione_timestep = leggi_pezzo_intestazione(0,t0); natoms=t0.natoms(); - buffer_tipi=new int [natoms]; - buffer_tipi_id=new int [natoms]; + buffer_type=new int [natoms]; + buffer_type_id=new int [natoms]; for (size_t i=0;i -size_t Traiettoria::leggi_pezzo(const size_t &partenza /// offset da cui partire (in "file") +size_t Trajectory::leggi_pezzo(const size_t &partenza /// offset da cui partire (in "file") ,TimestepManager ×tep /// timesteps header ,Chunk * &chunk /// actual data of the trajectory ){ @@ -281,11 +281,11 @@ size_t Traiettoria::leggi_pezzo(const size_t &partenza /// offset da cui partire if (timestep.triclinic()) { triclinic=true; - buffer_scatola_stride=9; + buffer_boxes_stride=9; box_format=BoxFormat::Lammps_triclinic; } else { triclinic=false; - buffer_scatola_stride=6; + buffer_boxes_stride=6; box_format=BoxFormat::Lammps_ortho; } @@ -321,7 +321,7 @@ size_t Traiettoria::leggi_pezzo(const size_t &partenza /// offset da cui partire * restituisce la dimensione in byte del pezzo letto, comprensiva dei chunk con i dati veri e propri che però non vengono restituiti * (quindi il pezzo successivo si troverà ad un offset di "partenza" + "il risultato di questa chiamata") **/ -size_t Traiettoria::leggi_pezzo_intestazione(const size_t &partenza /// offset da cui partire (in "file") +size_t Trajectory::leggi_pezzo_intestazione(const size_t &partenza /// offset da cui partire (in "file") ,TimestepManager ×tep /// qui viene restituito un oggetto TimestepManager ){ Chunk * ch; @@ -329,7 +329,7 @@ size_t Traiettoria::leggi_pezzo_intestazione(const size_t &partenza /// offset d } -Traiettoria::Errori Traiettoria::imposta_dimensione_finestra_accesso(const size_t &Ntimesteps){ +Trajectory::Errori Trajectory::set_data_access_block_size(const size_t &Ntimesteps){ if (!ok) { std::cerr << "mmap not correctly initialized!\n"; return non_inizializzato; @@ -342,19 +342,19 @@ Traiettoria::Errori Traiettoria::imposta_dimensione_finestra_accesso(const size_ if (loaded_timesteps!=Ntimesteps){ //questo alloca la memoria in modo corretto per permettere l'utilizzo delle istruzioni SIMD in fftw3 - fftw_free(buffer_posizioni); - fftw_free(buffer_posizioni_cm); - fftw_free(buffer_velocita_cm); - fftw_free(buffer_velocita); - buffer_posizioni_size=sizeof(double)*natoms*3*Ntimesteps; + fftw_free(buffer_positions); + fftw_free(buffer_positions_cm); + fftw_free(buffer_velocity_cm); + fftw_free(buffer_velocity); + buffer_positions_size=sizeof(double)*natoms*3*Ntimesteps; buffer_cm_size=sizeof(double)*3*Ntimesteps*ntypes; - buffer_posizioni= (double*) fftw_malloc(buffer_posizioni_size); - buffer_posizioni_cm= (double*) fftw_malloc(buffer_cm_size); - buffer_velocita_cm= (double*) fftw_malloc(buffer_cm_size); - buffer_velocita=(double*) fftw_malloc(buffer_posizioni_size); + buffer_positions= (double*) fftw_malloc(buffer_positions_size); + buffer_positions_cm= (double*) fftw_malloc(buffer_cm_size); + buffer_velocity_cm= (double*) fftw_malloc(buffer_cm_size); + buffer_velocity=(double*) fftw_malloc(buffer_positions_size); - delete [] buffer_scatola; - buffer_scatola= new double [buffer_scatola_stride*Ntimesteps]; + delete [] buffer_boxes; + buffer_boxes= new double [buffer_boxes_stride*Ntimesteps]; } loaded_timesteps=Ntimesteps; @@ -366,7 +366,7 @@ Traiettoria::Errori Traiettoria::imposta_dimensione_finestra_accesso(const size_ * alloca un array (timesteps) piu' lungo e copia il contenuto del vecchio indice in quello nuovo. * Inizializza a zero gli elementi in piu' che ci sono alla fine dell'array **/ -void Traiettoria::allunga_timesteps(size_t nuova_dimensione){ +void Trajectory::allunga_timesteps(size_t nuova_dimensione){ if (nuova_dimensione<=n_timesteps) return; size_t *tmp=new size_t [nuova_dimensione]; @@ -387,7 +387,7 @@ void Traiettoria::allunga_timesteps(size_t nuova_dimensione){ } -void Traiettoria::index_all() { +void Trajectory::index_all() { if (timestep_indicizzato>=n_timesteps-1) return; @@ -430,7 +430,7 @@ void Traiettoria::index_all() { #endif } -Traiettoria::Errori Traiettoria::imposta_inizio_accesso(const size_t ×tep) { +Trajectory::Errori Trajectory::set_access_at(const size_t ×tep) { cronometro cron; cron.start(); @@ -440,7 +440,7 @@ Traiettoria::Errori Traiettoria::imposta_inizio_accesso(const size_t ×tep) return non_inizializzato; } - if (timestep==timestep_corrente && dati_caricati) + if (timestep==current_timestep && dati_caricati) return Ok; #ifdef DEBUG @@ -453,7 +453,7 @@ Traiettoria::Errori Traiettoria::imposta_inizio_accesso(const size_t ×tep) //(vuol dire che la stima del numero di timesteps svolta all'inizio era sbagliata) if (timestep+loaded_timesteps>n_timesteps) { - if (fsize-timesteps[timestep_corrente]<(timestep-timestep_corrente + loaded_timesteps)*tstep_size) { + if (fsize-timesteps[current_timestep]<(timestep-current_timestep + loaded_timesteps)*tstep_size) { allunga_timesteps(timestep+loaded_timesteps+1); } else { std::cerr << "Cannot use timesteps beyond the end of file (access window was too big?)\n"; @@ -483,7 +483,7 @@ Traiettoria::Errori Traiettoria::imposta_inizio_accesso(const size_t ×tep) perror(nullptr); } - if (timestep < timestep_corrente && timesteps[timestep]+tstep_size*loaded_timesteps*2 < fsize ) { // avviso che le zone di memoria successive non mi servono + if (timestep < current_timestep && timesteps[timestep]+tstep_size*loaded_timesteps*2 < fsize ) { // avviso che le zone di memoria successive non mi servono res=madvise(file+timesteps[timestep]+tstep_size*loaded_timesteps*2, fsize-timesteps[timestep]+tstep_size*loaded_timesteps*2,MADV_DONTNEED); if ( res==-1) { @@ -525,25 +525,25 @@ Traiettoria::Errori Traiettoria::imposta_inizio_accesso(const size_t ×tep) size_t timestep_copy_tstart=0,timestep_copy_tend=0, timestep_read_start=timestep, timestep_read_end=timestep+loaded_timesteps; - ssize_t finestra_differenza=timestep_corrente-timestep; + ssize_t finestra_differenza=current_timestep-timestep; if (dati_caricati) { if (abs(finestra_differenza)0){ timestep_copy_tstart=0; timestep_copy_tend=loaded_timesteps-finestra_differenza; - timestep_read_end=timestep_corrente; + timestep_read_end=current_timestep; timestep_read_start=timestep; } else { timestep_copy_tstart=abs(finestra_differenza); timestep_copy_tend=loaded_timesteps; - timestep_read_start=timestep_corrente+loaded_timesteps; + timestep_read_start=current_timestep+loaded_timesteps; timestep_read_end=timestep+loaded_timesteps; } #ifdef DEBUG - std::cerr << "Copying data from buffer from "< #include "lammps_struct.h" #include "buffer_utils.h" #include "triclinic.h" -Traiettoria_numpy::Traiettoria_numpy(pybind11::buffer buffer_pos_, +Trajectory_numpy::Trajectory_numpy(pybind11::buffer buffer_pos_, pybind11::buffer buffer_vel_, pybind11::buffer buffer_types_, pybind11::buffer buffer_box_, - TraiettoriaBase::BoxFormat matrix_box, + BaseTrajectory::BoxFormat matrix_box, bool wrap_pbc_, bool save_rotation_matrix) : buffer_pos{buffer_pos_},buffer_vel{buffer_vel_},buffer_types{buffer_types_},buffer_box{buffer_box_},rotation_matrix{nullptr} { @@ -83,20 +83,20 @@ Traiettoria_numpy::Traiettoria_numpy(pybind11::buffer buffer_pos_, natoms=info_pos.shape[1]; n_timesteps=info_pos.shape[0]; - buffer_scatola_stride = 6; + buffer_boxes_stride = 6; triclinic=false; std::vector >> cells_qr; // first step, rotation matrix if (matrix_box== BoxFormat::Lammps_ortho || matrix_box==BoxFormat::Lammps_triclinic){ //convert the format in the internal one, a little more convenient for min image algorithm if (matrix_box==BoxFormat::Lammps_triclinic) { - buffer_scatola_stride = 9; + buffer_boxes_stride = 9; triclinic=true; } - buffer_scatola = new double[info_box.shape[0]*buffer_scatola_stride]; + buffer_boxes = new double[info_box.shape[0]*buffer_boxes_stride]; box_allocated=true; //copy and do a permutation of everything for (ssize_t i=0;i(info_box.ptr)[i*buffer_scatola_stride],buffer_scatola_stride*sizeof (double)); - lammps_to_internal(buffer_scatola+i*buffer_scatola_stride); + std::memcpy(buffer_boxes+i*buffer_boxes_stride,&static_cast(info_box.ptr)[i*buffer_boxes_stride],buffer_boxes_stride*sizeof (double)); + lammps_to_internal(buffer_boxes+i*buffer_boxes_stride); } std::cerr << "Input format is lammps"<(info_vel.ptr); + buffer_velocity=static_cast(info_vel.ptr); box_format=BoxFormat::Lammps_ortho; } else { - buffer_scatola_stride = 9; + buffer_boxes_stride = 9; box_format=BoxFormat::Lammps_triclinic; } if (wrap_pbc || triclinic){ // when I need also a copy of the positions - buffer_posizioni=new double[info_pos.shape[0]*info_pos.shape[1]*info_pos.shape[2]]; + buffer_positions=new double[info_pos.shape[0]*info_pos.shape[1]*info_pos.shape[2]]; posizioni_allocated=true; }else { posizioni_allocated=false; - buffer_posizioni=static_cast(info_pos.ptr); + buffer_positions=static_cast(info_pos.ptr); velocities_allocated=false; - buffer_velocita=static_cast(info_vel.ptr); + buffer_velocity=static_cast(info_vel.ptr); } //now everything is allocated/moved. Do the work of translation to the lammps (wapped) format @@ -161,8 +161,8 @@ Traiettoria_numpy::Traiettoria_numpy(pybind11::buffer buffer_pos_, for (ssize_t i=0;i(info_pos.ptr)[i*3*natoms+n*3],3*sizeof (double)); std::memcpy(vel,&static_cast(info_vel.ptr)[i*3*natoms+n*3],3*sizeof (double)); cells_qr[cur_idx].second.rotate_vec(pos); @@ -171,16 +171,16 @@ Traiettoria_numpy::Traiettoria_numpy(pybind11::buffer buffer_pos_, if (cur_idx + 1 < cells_qr.size() && cells_qr[cur_idx+1].first-1 == i) cur_idx++; } } else if (wrap_pbc) { - std::memcpy(buffer_posizioni,info_pos.ptr,sizeof(double)*3*natoms*n_timesteps); + std::memcpy(buffer_positions,info_pos.ptr,sizeof(double)*3*natoms*n_timesteps); } - buffer_tipi=static_cast(info_types.ptr); - buffer_tipi_id=new int[natoms]; + buffer_type=static_cast(info_types.ptr); + buffer_type_id=new int[natoms]; get_ntypes(); //init type ids // //calculate center of mass velocity and position (without pbc) - calc_cm_pos_vel(static_cast(info_pos.ptr),buffer_posizioni_cm); - calc_cm_pos_vel(static_cast(info_vel.ptr),buffer_velocita_cm); + calc_cm_pos_vel(static_cast(info_pos.ptr),buffer_positions_cm); + calc_cm_pos_vel(static_cast(info_vel.ptr),buffer_velocity_cm); if (wrap_pbc) { @@ -199,7 +199,7 @@ Traiettoria_numpy::Traiettoria_numpy(pybind11::buffer buffer_pos_, } void -Traiettoria_numpy::calc_cm_pos_vel(double * a, double * & cm){ +Trajectory_numpy::calc_cm_pos_vel(double * a, double * & cm){ if (cm==nullptr){ cm=new double[3*n_timesteps*ntypes]; } else { @@ -224,17 +224,17 @@ Traiettoria_numpy::calc_cm_pos_vel(double * a, double * & cm){ delete [] cont; } -Traiettoria_numpy::~Traiettoria_numpy() { +Trajectory_numpy::~Trajectory_numpy() { if (box_allocated) - delete [] buffer_scatola; + delete [] buffer_boxes; if (velocities_allocated) - delete [] buffer_velocita; + delete [] buffer_velocity; if (posizioni_allocated) - delete [] buffer_posizioni; - delete [] buffer_tipi_id; - delete [] masse; - delete [] cariche; - delete [] buffer_posizioni_cm; - delete [] buffer_velocita_cm; + delete [] buffer_positions; + delete [] buffer_type_id; + delete [] mass; + delete [] charge; + delete [] buffer_positions_cm; + delete [] buffer_velocity_cm; delete [] rotation_matrix; } diff --git a/pyanalisi/src/pyanalisi.cpp b/pyanalisi/src/pyanalisi.cpp index eb3b973db..a846f2873 100644 --- a/pyanalisi/src/pyanalisi.cpp +++ b/pyanalisi/src/pyanalisi.cpp @@ -3,17 +3,14 @@ #include "pybind11/pybind11.h" #include "pybind11/stl.h" #include -#include "traiettoria.h" +#include "trajectory.h" #include "spettrovibrazionale.h" #include "readlog.h" #include "correlatorespaziale.h" #include "gofrt.h" #include "greenkuboNcomponentionicfluid.h" -#include "heatc.h" -#include "centerdiff.h" -#include "centerofmassdiff.h" #include "msd.h" -#include "traiettoria_numpy.h" +#include "trajectory_numpy.h" #include "readlog_numpy.h" #include "sphericalcorrelations.h" #include "config.h" @@ -43,18 +40,18 @@ provide a list of [(max_neighbours, rmax^2, rmax_2^2),...] if you want to use th if true compute the averaged steinhardt order parameter (an additional loop over neighbours is performed). Implemented only with SANN )lol") .def("reset",&SOP::reset) - .def("calculate", &SOP::calcola) + .def("calculate", &SOP::calculate) .def_buffer([](SOP & m) -> py::buffer_info { /* std::cerr <<"shape ("<< m.get_shape().size() << "): "; for (auto & n: m.get_shape()) std::cerr << n << " "; std::cerr <::format(), m.get_shape().size(), @@ -72,11 +69,11 @@ void define_atomic_traj(py::module & m, std::string typestr){ calculates g(r) and, in general, g(r,t). Interesting dynamical property Parameters ---------- Trajectory instance rmin rmax nbin maximum time lag number of threads time skip every time debug flag )begend") .def("reset",& Gofrt::reset,R"begend()begend") - .def("getNumberOfExtraTimestepsNeeded",&Gofrt::numeroTimestepsOltreFineBlocco ,R"begend()begend") - .def("calculate", & Gofrt::calcola,R"begend()begend") + .def("getNumberOfExtraTimestepsNeeded",&Gofrt::nExtraTimesteps ,R"begend()begend") + .def("calculate", & Gofrt::calculate,R"begend()begend") .def_buffer([](Gofrt & g) -> py::buffer_info { return py::buffer_info( - g.accesso_lista(), /* Pointer to buffer */ + g.access_vdata(), /* Pointer to buffer */ sizeof(double), /* Size of one scalar */ py::format_descriptor::format(), /* Python struct-style format descriptor */ g.get_shape().size(), /* Number of dimensions */ @@ -101,10 +98,10 @@ void define_atomic_traj(py::module & m, std::string typestr){ provide a list of [(max_neighbours, rmax^2, rmax_2^2),...] if you want to use the SANN algorithm and use ibin=0 always )lol") .def("reset",&SHC::reset) - .def("calculate", &SHC::calcola) + .def("calculate", &SHC::calculate) .def_buffer([](SHC & m) -> py::buffer_info { return py::buffer_info( - m.accesso_lista(), + m.access_vdata(), sizeof(double), py::format_descriptor::format(), m.get_shape().size(), @@ -133,11 +130,11 @@ void define_atomic_traj(py::module & m, std::string typestr){ debug flag )begend") .def("reset",&MSD::reset) - .def("getNumberOfExtraTimestepsNeeded",&MSD::numeroTimestepsOltreFineBlocco) - .def("calculate", &MSD::calcola) + .def("getNumberOfExtraTimestepsNeeded",&MSD::nExtraTimesteps) + .def("calculate", &MSD::calculate) .def_buffer([](MSD & m) -> py::buffer_info { return py::buffer_info( - m.accesso_lista(), + m.access_vdata(), sizeof(double), py::format_descriptor::format(), m.get_shape().size(), @@ -152,10 +149,10 @@ void define_atomic_traj(py::module & m, std::string typestr){ py::class_(m,(std::string("PositionHistogram")+typestr).c_str(),py::buffer_protocol()) .def(py::init,unsigned, unsigned>()) .def("reset",&AD::reset) - .def("calculate",&AD::calcola) + .def("calculate",&AD::calculate) .def_buffer([](AD & m)->py::buffer_info { return py::buffer_info( - m.accesso_lista(), + m.access_vdata(), sizeof(long), py::format_descriptor::format(), m.get_shape().size(), @@ -174,11 +171,11 @@ void define_atomic_traj(py::module & m, std::string typestr){ debug flag )begend") .def("reset",&VDOS::reset) - .def("calculate", &VDOS::calcola) + .def("calculate", &VDOS::calculate) .def("spectrum", &VDOS::spettro) .def_buffer([](VDOS & m) -> py::buffer_info { return py::buffer_info( - m.accesso_lista(), + m.access_vdata(), sizeof(double), py::format_descriptor::format(), m.get_shape().size(), @@ -258,14 +255,14 @@ void gk(py::module & m, std::string typestr){ unsigned int -> benchmark parameter start unsigned int -> benchmark parameter stop )begend") - .def("getNumberOfExtraTimestepsNeeded",&GK::numeroTimestepsOltreFineBlocco,R"begend( + .def("getNumberOfExtraTimestepsNeeded",&GK::nExtraTimesteps,R"begend( perform an efficient, multithreaded Green-Kubo calculations of the transport coefficient of a condensed matter system, given the currents. )begend") .def("reset",&GK::reset,R"begend()begend") - .def("calculate",&GK::calcola,R"begend()begend") + .def("calculate",&GK::calculate,R"begend()begend") .def_buffer([](GK & g) -> py::buffer_info { return py::buffer_info( - g.accesso_lista(), /* Pointer to buffer */ + g.access_vdata(), /* Pointer to buffer */ sizeof(double), /* Size of one scalar */ py::format_descriptor::format(), /* Python struct-style format descriptor */ g.get_shape().size(), /* Number of dimensions */ @@ -287,13 +284,13 @@ R"lol( )lol") .def("get_positions_copy", [](Tk & t) { double * foo=nullptr; - if (t.posizioni(t.get_current_timestep(),0) == nullptr) { + if (t.positions(t.get_current_timestep(),0) == nullptr) { return pybind11::array_t(); } long nts=t.get_nloaded_timesteps(); long nat=t.get_natoms(); foo = new double[nts*nat*3]; - std::memcpy(foo,t.posizioni(t.get_current_timestep(),0),sizeof (double)*nts*nat*3); + std::memcpy(foo,t.positions(t.get_current_timestep(),0),sizeof (double)*nts*nat*3); pybind11::capsule free_when_done(foo, [](void *f) { double *foo = reinterpret_cast(f); std::cerr << "freeing memory @ " << f << "\n"; @@ -306,14 +303,14 @@ R"lol( free_when_done );}) .def("get_velocities_copy", [](Tk & t) { - double * foo=t.velocita(t.get_current_timestep(),0); + double * foo=t.velocity(t.get_current_timestep(),0); if (foo == nullptr) { return pybind11::array_t(); } long nts=t.get_nloaded_timesteps(); long nat=t.get_natoms(); foo = new double[nts*nat*3]; - std::memcpy(foo,t.velocita(t.get_current_timestep(),0),sizeof (double)*nts*nat*3); + std::memcpy(foo,t.velocity(t.get_current_timestep(),0),sizeof (double)*nts*nat*3); pybind11::capsule free_when_done(foo, [](void *f) { double *foo = reinterpret_cast(f); std::cerr << "freeing memory @ " << f << "\n"; @@ -326,14 +323,14 @@ R"lol( free_when_done );}) .def("get_box_copy", [](Tk & t) { - double * foo=t.scatola(t.get_current_timestep()); + double * foo=t.box(t.get_current_timestep()); if (foo == nullptr) { return pybind11::array_t(); } long nts=t.get_nloaded_timesteps(); long nb=t.get_box_stride(); foo = new double[nts*nb]; - std::memcpy(foo,t.scatola(t.get_current_timestep()),sizeof (double)*nts*nb); + std::memcpy(foo,t.box(t.get_current_timestep()),sizeof (double)*nts*nb); pybind11::capsule free_when_done(foo, [](void *f) { double *foo = reinterpret_cast(f); std::cerr << "freeing memory @ " << f << "\n"; @@ -370,21 +367,21 @@ R"lol( PYBIND11_MODULE(pyanalisi,m) { #ifdef BUILD_MMAP - trajectory_common_interfaces>( - py::class_(m,"Traj", py::buffer_protocol())) + trajectory_common_interfaces>( + py::class_(m,"Traj", py::buffer_protocol())) .def(py::init(),R"begend( Parameters ---------- name of the binary file to open )begend") - .def("setWrapPbc",&Traiettoria::set_pbc_wrap,R"begend( + .def("setWrapPbc",&Trajectory::set_pbc_wrap,R"begend( wrap all the atomic coordinate inside the simulation box, red from the binary file Parameters ---------- True/False )begend") - .def("setAccessWindowSize",[](Traiettoria & t,int ts) {return (int) t.imposta_dimensione_finestra_accesso(ts);},R"begend( + .def("setAccessWindowSize",[](Trajectory & t,int ts) {return (int) t.set_data_access_block_size(ts);},R"begend( sets the size of the read block. Must fit in memory Parameters @@ -396,27 +393,27 @@ PYBIND11_MODULE(pyanalisi,m) { 1 -> success 0/-1 -> failure )begend") - .def("setAccessStart",[](Traiettoria & t,int ts) {return (int) t.imposta_inizio_accesso(ts);},R"begend( + .def("setAccessStart",[](Trajectory & t,int ts) {return (int) t.set_access_at(ts);},R"begend( sets the first timestep to read, and read the full block Parameters ---------- int -> first timestep to read )begend") - //.def("getPosition",[](Traiettoria & t, int ts,int natoms){return std::array<>}) - .def("serving_pos", &Traiettoria::serving_pos) - .def("toggle_pos_vel", &Traiettoria::toggle_pos_vel) - .def_buffer([](Traiettoria & g) -> py::buffer_info { + //.def("getPosition",[](Trajectory & t, int ts,int natoms){return std::array<>}) + .def("serving_pos", &Trajectory::serving_pos) + .def("toggle_pos_vel", &Trajectory::toggle_pos_vel) + .def_buffer([](Trajectory & g) -> py::buffer_info { return py::buffer_info( - g.serving_pos()? g.posizioni_inizio() : g.velocita_inizio(), /* Pointer to buffer */ + g.serving_pos()? g.positions_data() : g.velocity_data(), /* Pointer to buffer */ sizeof(double), /* Size of one scalar */ py::format_descriptor::format(), /* Python struct-style format descriptor */ g.get_shape().size(), /* Number of dimensions */ g.get_shape(), /* Buffer dimensions */ g.get_stride()); }) - .def("get_lammps_id", [](Traiettoria & t) { + .def("get_lammps_id", [](Trajectory & t) { int * foo=t.get_lammps_id(); pybind11::capsule free_when_done(foo, [](void *f) { int *foo = reinterpret_cast(f); @@ -430,7 +427,7 @@ PYBIND11_MODULE(pyanalisi,m) { foo, free_when_done );}) - .def("get_lammps_type", [](Traiettoria & t) { + .def("get_lammps_type", [](Trajectory & t) { int * foo=t.get_lammps_type(); pybind11::capsule free_when_done(foo, [](void *f) { int *foo = reinterpret_cast(f); @@ -450,7 +447,7 @@ PYBIND11_MODULE(pyanalisi,m) { using RL = ReadLog; py::class_(m,"ReadLog_file") - .def(py::init >() ,R"begend( + .def(py::init >() ,R"begend( perform an advanced read of a column formatted textfile with headers. Support multithreaded reading of chunks of lines. Parameters @@ -480,9 +477,9 @@ PYBIND11_MODULE(pyanalisi,m) { return number of timesteps read )begend"); - trajectory_common_interfaces>( - py::class_(m,"Trajectory")) - .def(py::init(), + trajectory_common_interfaces>( + py::class_(m,"Trajectory")) + .def(py::init(), py::keep_alive<1, 2>(), py::keep_alive<1, 3>(), py::keep_alive<1, 4>(), @@ -501,7 +498,7 @@ PYBIND11_MODULE(pyanalisi,m) { - .def("get_rotation_matrix", [](Traiettoria_numpy & t) { + .def("get_rotation_matrix", [](Trajectory_numpy & t) { double * foo=t.get_rotation_matrix(0); if (foo == nullptr) { return pybind11::array_t(); @@ -521,17 +518,17 @@ PYBIND11_MODULE(pyanalisi,m) { free_when_done );}) ; - py::enum_(m, "BoxFormat", py::arithmetic()) - .value("Invalid", Traiettoria_numpy::BoxFormat::Invalid) - .value("CellVectors", Traiettoria_numpy::BoxFormat::Cell_vectors) - .value("LammpsOrtho", Traiettoria_numpy::BoxFormat::Lammps_ortho) - .value("LammpsTriclinic", Traiettoria_numpy::BoxFormat::Lammps_triclinic); + py::enum_(m, "BoxFormat", py::arithmetic()) + .value("Invalid", Trajectory_numpy::BoxFormat::Invalid) + .value("CellVectors", Trajectory_numpy::BoxFormat::Cell_vectors) + .value("LammpsOrtho", Trajectory_numpy::BoxFormat::Lammps_ortho) + .value("LammpsTriclinic", Trajectory_numpy::BoxFormat::Lammps_triclinic); #ifdef BUILD_MMAP - define_atomic_traj(m,"_lammps"); + define_atomic_traj(m,"_lammps"); gk >(m,"_file"); #endif - define_atomic_traj(m,""); + define_atomic_traj(m,""); gk >(m,""); m.def("info",[]() -> std::string{ return _info_msg ; diff --git a/tests/src/test_fixtures.h b/tests/src/test_fixtures.h index b21783197..ef2e87c6b 100644 --- a/tests/src/test_fixtures.h +++ b/tests/src/test_fixtures.h @@ -5,7 +5,7 @@ #include #include #include "config.h" -#include "traiettoria.h" +#include "trajectory.h" #include "readlog.h" struct TestPath{ @@ -16,11 +16,11 @@ struct TestPath{ struct TrajSetup{ explicit TrajSetup(bool pbc=false,bool f2020=false): traj{f2020? path.path+"data/lammps2020.bin" : path.path+"data/lammps.bin"}{ traj.set_pbc_wrap(pbc); - traj.imposta_dimensione_finestra_accesso(150); - traj.imposta_inizio_accesso(0); + traj.set_data_access_block_size(150); + traj.set_access_at(0); } const TestPath path; - Traiettoria traj; + Trajectory traj; }; diff --git a/tests/src/test_gk.cpp b/tests/src/test_gk.cpp index 42d81734c..1c5572a36 100644 --- a/tests/src/test_gk.cpp +++ b/tests/src/test_gk.cpp @@ -20,16 +20,16 @@ struct GkFixture{ BOOST_FIXTURE_TEST_SUITE(gk_traj_1, GkFixture<1>) BOOST_AUTO_TEST_CASE(test_gk_2comp_nth1) { - gk.calcola(0); - BOOST_TEST(data.test_regression("gk_2comp",gk.accesso_lista(),gk.lunghezza())); + gk.calculate(0); + BOOST_TEST(data.test_regression("gk_2comp",gk.access_vdata(),gk.lunghezza())); } BOOST_AUTO_TEST_SUITE_END() BOOST_FIXTURE_TEST_SUITE(gk_traj_2, GkFixture<2>) BOOST_AUTO_TEST_CASE(test_gk_2comp_nth2) { - gk.calcola(0); - BOOST_TEST(data.test_regression("gk_2comp",gk.accesso_lista(),gk.lunghezza())); + gk.calculate(0); + BOOST_TEST(data.test_regression("gk_2comp",gk.access_vdata(),gk.lunghezza())); } BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/src/test_lammps2020.cpp b/tests/src/test_lammps2020.cpp index c5703347c..646e867e2 100644 --- a/tests/src/test_lammps2020.cpp +++ b/tests/src/test_lammps2020.cpp @@ -11,11 +11,11 @@ struct MsdFixture { msd{&traj.traj, 1,0,NTH,true} {} TrajSetup traj; - MSD msd; + MSD msd; DataRegression data; - MSD & calc(int primo) { + MSD & calc(int primo) { msd.reset(75-primo); - msd.calcola(primo); + msd.calculate(primo); return msd; } size_t size(){auto s= msd.get_shape(); return s[0]*s[1]*s[2];} @@ -37,11 +37,11 @@ struct GofrFixture { } {} TrajSetup traj; - Gofrt gofr; + Gofrt gofr; DataRegression data; - Gofrt & calc(int primo) { + Gofrt & calc(int primo) { gofr.reset(75-primo); - gofr.calcola(primo); + gofr.calculate(primo); return gofr; } size_t size(){auto s= gofr.get_shape(); return s[0]*s[1]*s[2];} @@ -50,8 +50,8 @@ struct GofrFixture { BOOST_AUTO_TEST_CASE(test_cell_permutator){ double d[6]={1,2,3,4,5,6}; double e[6]={1,2,3,4,5,6}; - Traiettoria::lammps_to_internal(d); - Traiettoria::internal_to_lammps(d); + Trajectory::lammps_to_internal(d); + Trajectory::internal_to_lammps(d); for (int i=0;i<6;++i){ BOOST_TEST(d[i]==e[i]); } @@ -62,9 +62,9 @@ using MsdFixture_ ## N = MsdFixture;\ BOOST_FIXTURE_TEST_SUITE(test_msd##N, MsdFixture_ ## N)\ BOOST_AUTO_TEST_CASE(test_msd_##N){\ calc(0);\ - BOOST_TEST(data.test_regression("msd_"#N"a",msd.accesso_lista(),size()));\ + BOOST_TEST(data.test_regression("msd_"#N"a",msd.access_vdata(),size()));\ calc(13);\ - BOOST_TEST(data.test_regression("msd_"#N"b",msd.accesso_lista(),size()));\ + BOOST_TEST(data.test_regression("msd_"#N"b",msd.access_vdata(),size()));\ }\ BOOST_AUTO_TEST_SUITE_END() @@ -78,9 +78,9 @@ using GofrFixture_ ## N ## _ ## T = GofrFixture;\ BOOST_FIXTURE_TEST_SUITE(test_gofr##N ## _ ## T, GofrFixture_ ## N ## _ ## T)\ BOOST_AUTO_TEST_CASE(test_gofr_##N ## _ ## T){\ calc(0);\ - BOOST_TEST(data.test_regression("gofr_"#N"a"#T,gofr.accesso_lista(),size()));\ + BOOST_TEST(data.test_regression("gofr_"#N"a"#T,gofr.access_vdata(),size()));\ calc(13);\ - BOOST_TEST(data.test_regression("gofr_"#N"b"#T,gofr.accesso_lista(),size()));\ + BOOST_TEST(data.test_regression("gofr_"#N"b"#T,gofr.access_vdata(),size()));\ }\ BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/src/test_position_histogram.cpp b/tests/src/test_position_histogram.cpp index baf9747e1..26254c518 100644 --- a/tests/src/test_position_histogram.cpp +++ b/tests/src/test_position_histogram.cpp @@ -10,13 +10,13 @@ struct PhFixture { data.path+="ad/"; } TrajSetup traj; - using AD = AtomicDensity; + using AD = AtomicDensity; AD ret; DataRegression data; AD & calc(int primo){ ret.reset(150-primo); - ret.calcola(primo); + ret.calculate(primo); return ret; } size_t size(){auto s= ret.get_shape(); return s[0]*s[1]*s[2];} @@ -28,9 +28,9 @@ BOOST_FIXTURE_TEST_SUITE(test_atomic_density ## N, PhFixture_ ## N )\ BOOST_AUTO_TEST_CASE(test_position_histogram ## N)\ {\ calc(0);\ - BOOST_TEST(data.test_regression("test_calcola_0",ret.accesso_lista(),size()));\ + BOOST_TEST(data.test_regression("test_calcola_0",ret.access_vdata(),size()));\ calc(12);\ - BOOST_TEST(data.test_regression("test_calcola_1",ret.accesso_lista(),size()));\ + BOOST_TEST(data.test_regression("test_calcola_1",ret.access_vdata(),size()));\ }\ BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/src/test_sh.cpp b/tests/src/test_sh.cpp index 36229d1d2..f3c224d62 100644 --- a/tests/src/test_sh.cpp +++ b/tests/src/test_sh.cpp @@ -18,7 +18,7 @@ struct ShFixture{ const unsigned int nbin; const size_t natoms; const size_t ntypes; - SphericalCorrelations sh; + SphericalCorrelations sh; DataRegression data; double workspace[(l+1)*(l+1)], cheby[(l+1)*2]; void calc(int timestep, double * res){ @@ -43,7 +43,7 @@ struct SteinhardtFixture{ const unsigned int nbin; const size_t natoms; const size_t ntypes; - Steinhardt sh; + Steinhardt sh; DataRegression data; }; template @@ -57,7 +57,7 @@ struct SteinhardtFixtureNeigh{ const unsigned int nbin; const size_t natoms; const size_t ntypes; - Steinhardt sh; + Steinhardt sh; DataRegression data; }; @@ -80,8 +80,8 @@ typedef SteinhardtFixtureNeigh<6,3> StFix6_tr_3; BOOST_FIXTURE_TEST_SUITE(sh ## T, T )\ BOOST_AUTO_TEST_CASE(test_calcola)\ {\ - sh.calcola(0);\ - BOOST_TEST(data.test_regression("test_calcola" # suff ,sh.accesso_lista(),sh.lunghezza()));\ + sh.calculate(0);\ + BOOST_TEST(data.test_regression("test_calcola" # suff ,sh.access_vdata(),sh.lunghezza()));\ }\ BOOST_AUTO_TEST_SUITE_END() diff --git a/tests/src/test_trajectory.cpp b/tests/src/test_trajectory.cpp index 3d95f85e9..0451b2038 100644 --- a/tests/src/test_trajectory.cpp +++ b/tests/src/test_trajectory.cpp @@ -14,7 +14,7 @@ struct NeighTest{ NeighTest(bool l2020=true):t{true,l2020},n{&t.traj,{{99,1.0,1.0},{67,1.0,1.0}}} {}; TrajSetup t; DataRegression data; - Neighbours n; + Neighbours n; }; @@ -28,7 +28,7 @@ BOOST_AUTO_TEST_CASE(min_image){ for (size_t j=0;j