diff --git a/AGENTS.md b/AGENTS.md index ea0fcd28fc7..45600683201 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -89,6 +89,16 @@ rules. Read the complete governance document before making or reviewing changes: `docs/advanced/`, `toolchain/`, `Dockerfile.gnu`, `Dockerfile.intel`, and `Dockerfile.cuda`. +## File Move / Split Checklist + +When moving, splitting, or deleting a `.cpp`/`.h`, update **all** build surfaces: + +1. `source/Makefile.Objects` (legacy Makefile object lists) +2. All C++ `CMakeLists.txt` +3. `python/pyabacus/src/*/CMakeLists.txt` — pyabacus has independent hardcoded source lists that do not track the main CMake build + +When splitting a function into a new `.cpp`, if its callers span LCAO and PW, put it in an unconditionally-built OBJECT library, not inside `if(ENABLE_LCAO)`. + ## Build And Test Entry Points - Prefer the repository CMake/CTest flow already used by CI. For focused local diff --git a/python/pyabacus/src/ModuleNAO/CMakeLists.txt b/python/pyabacus/src/ModuleNAO/CMakeLists.txt index b8ee95319cf..0366d1c9122 100644 --- a/python/pyabacus/src/ModuleNAO/CMakeLists.txt +++ b/python/pyabacus/src/ModuleNAO/CMakeLists.txt @@ -8,6 +8,7 @@ list(APPEND _naos ${NAO_PATH}/radial_collection.cpp ${NAO_PATH}/radial_set.cpp ${NAO_PATH}/real_gaunt_table.cpp + ${NAO_PATH}/smoothgen.cpp ${NAO_PATH}/sphbes_radials.cpp ${NAO_PATH}/two_center_bundle.cpp ${NAO_PATH}/two_center_integrator.cpp diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 05ee4118ee5..2322dda4976 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -733,6 +733,7 @@ target_link_libraries( elecstate module_pwdft module_proj + nao_smoothgen module_ofdft module_stodft module_dfpt diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 7514f75060a..839b87d20bb 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -144,7 +144,6 @@ OBJS_BASE=assoc_laguerre.o\ complexmatrix.o\ clebsch_gordan_coeff.o\ element_basis_index.o\ - export.o\ inverse_matrix.o\ gather_math_lib_info.o\ global_file.o\ @@ -557,7 +556,6 @@ OBJS_ORBITAL=orb_atomic.o\ two_center_bundle.o\ two_center_integrator.o\ two_center_table.o\ - projgen.o\ OBJS_PSI=psi.o\ @@ -992,6 +990,7 @@ OBJS_SRCPW=h_ewald_pw.o\ onsite_proj_tools_force.o\ onsite_proj_tools_stress.o\ onsite_proj_print.o\ + smoothgen.o\ vsep_pw.o OBJS_VDW=vdw.o\ diff --git a/source/source_base/CMakeLists.txt b/source/source_base/CMakeLists.txt index 4d305ad41d3..ccf5aae9077 100644 --- a/source/source_base/CMakeLists.txt +++ b/source/source_base/CMakeLists.txt @@ -17,7 +17,6 @@ add_library( complexarray.cpp complexmatrix.cpp element_basis_index.cpp - export.cpp inverse_matrix.cpp gather_math_lib_info.cpp global_file.cpp @@ -61,7 +60,6 @@ add_library( sph_bessel_tf.cpp cubic_spline.cpp parallel_2d.cpp - projgen.cpp module_out/binstream.cpp module_out/filename.cpp module_out/sparse_matrix.cpp diff --git a/source/source_base/complexmatrix.h b/source/source_base/complexmatrix.h index 9dc890198d3..c47ab8f333e 100644 --- a/source/source_base/complexmatrix.h +++ b/source/source_base/complexmatrix.h @@ -5,8 +5,6 @@ #include "matrix.h" -#ifdef _MCD_CHECK -#endif namespace ModuleBase { class ComplexMatrix diff --git a/source/source_base/export.cpp b/source/source_base/export.cpp deleted file mode 100644 index 947f12ffb91..00000000000 --- a/source/source_base/export.cpp +++ /dev/null @@ -1,33 +0,0 @@ -// ============================================================================= -// C++ Header File -// Project: -// File: export.h -// Author: mohan -// Comment: -// Warning: -// Start time: 2008-9-3 -// Last modified: -// ============================================================================= - -#include "export.h" -/* -void ModuleBase::IF_MATCH(const std::string &name,const std::string &name2) -{ - if(name!=name2) - { - if(GlobalV::MY_RANK == 0) - { - std::cout<<"\n Can not match : "< -#include -#include -#include -#include -#include "global_variable.h" - -#ifdef __MPI -#include "mpi.h" -#endif -//these two function is not used yet!!!! -/*namespace ModuleBase -{ - -template -void ModuleBase::GlobalFunc::AUTO_SET(std::ofstream &ofs,const std::string &name,const T &a) -{ - ofs<<" AUTO_SET "< -void IF_MATCH(const T &a,const T &b) -{ - if(a!=b) - { - if(GlobalV::MY_RANK == 0) - { - std::cout<<"\n Can not match : "< -#include -#include -#include - - -void calculate(); - -int main(int argc, char **argv) -{ - - std::cout << "Hello, this is the 'base' module of ABACUS." << std::endl; - - std::cout << "The module searchs for the neighboring atoms for a given atomic position" << std::endl; - - std::cout << "Right now, the module is still empty, soon we will have more tests." << std::endl; - - calculate(); - - return 0; -} - - -void calculate() -{ - std::ofstream ofs("log.txt"); - -// ooo.set_orb_tables(); - - ofs.close(); - - std::cout << "--------------------" << std::endl; - std::cout << " Have a great day! " << std::endl; - std::cout << "--------------------" << std::endl; - -/* - time_t time_start = std::time(NULL); - -// ModuleBase::timer::start(); - - //---------------------------------------------------------- - // main program for doing electronic structure calculations - //---------------------------------------------------------- -// Driver DD; -// DD.init(); - - time_t time_finish= std::time(NULL); - - // print out information before ABACUS ends - std::cout << "\n START Time : " << ctime(&time_start); - std::cout << " FINISH Time : " << ctime(&time_finish); - std::cout << " TOTAL Time : " << difftime(time_finish, time_start) << std::endl; - - double total_time = difftime(time_finish, time_start); - int hour = total_time / 3600; - int mins = ( total_time - 3600 * hour ) / 60; - int secs = total_time - 3600 * hour - 60 * mins ; - std::cout << " Total Time : " << hour << " h " - << mins << " mins " - << secs << " secs "<< std::endl; -*/ - - return; -} diff --git a/source/source_base/matrix.h b/source/source_base/matrix.h index 8d2df96c45b..866bd6164c1 100644 --- a/source/source_base/matrix.h +++ b/source/source_base/matrix.h @@ -3,9 +3,6 @@ // Peize Lin update 2018-07-29 -#ifdef _MCD_CHECK -#endif - #include #include diff --git a/source/source_base/matrix3.h b/source/source_base/matrix3.h index cfd08972716..953953fbd5a 100644 --- a/source/source_base/matrix3.h +++ b/source/source_base/matrix3.h @@ -1,10 +1,6 @@ #ifndef MatriX3_H #define MatriX3_H -#ifdef _MCD_CHECK -#include "mcd.h" -#endif - #include "matrix.h" #include "vector3.h" diff --git a/source/source_base/matrix_wrapper.h b/source/source_base/matrix_wrapper.h index 9745c8eb5f2..67f2d520d7e 100644 --- a/source/source_base/matrix_wrapper.h +++ b/source/source_base/matrix_wrapper.h @@ -3,8 +3,6 @@ // DATE : 2018-07-31 //========================================================== -#ifndef MATRIX_WRAPPER_TIANHE2 - #ifndef MATRIX_WRAPPER_H #define MATRIX_WRAPPER_H @@ -86,7 +84,3 @@ inline Matrix_Wrapper& Matrix_Wrapper::operator=( Matrix_Wrapper&&m ) } #endif - -#else -#include "matrix_wrapper_tianhe2.h" -#endif diff --git a/source/source_base/matrix_wrapper_tianhe2.h b/source/source_base/matrix_wrapper_tianhe2.h deleted file mode 100644 index 79d7adbb0ff..00000000000 --- a/source/source_base/matrix_wrapper_tianhe2.h +++ /dev/null @@ -1,55 +0,0 @@ -//========================================================== -// AUTHOR : Peize Lin -// DATE : 2018-07-31 -//========================================================== - -#ifndef MATRIX_WRAPPER_H -#define MATRIX_WRAPPER_H - -#include "matrix.h" -#include -namespace ModuleBase -{ - -// !!! Attention: c is very dangerous, may point to somewhere deleted - -class Matrix_Wrapper -{ -public: - double * c = nullptr; - int nr; - int nc; - - Matrix_Wrapper(): nr(0), nc(0), c(nullptr){} - Matrix_Wrapper( const matrix &m ): nr(m.nr), nc(m.nc), c(m.c){} - inline void create( const int nr_in, const int nc_in, const bool flag_zero ); - inline matrix to_matrix(); - Matrix_Wrapper( const Matrix_Wrapper &m ): nr(m.nr), nc(m.nc), c(m.c){} - Matrix_Wrapper( Matrix_Wrapper &m ): nr(m.nr), nc(m.nc), c(m.c){} - Matrix_Wrapper( Matrix_Wrapper &&m ): nr(m.nr), nc(m.nc), c(m.c){} - inline Matrix_Wrapper&operator=( const Matrix_Wrapper&m ){ nr=m.nr; nc=m.nc; c=m.c; return *this; }; - inline Matrix_Wrapper&operator=( Matrix_Wrapper&m ){ nr=m.nr; nc=m.nc; c=m.c; return *this; }; - inline Matrix_Wrapper&operator=( Matrix_Wrapper&&m ){ nr=m.nr; nc=m.nc; c=m.c; return *this; }; -}; - - -inline void Matrix_Wrapper::create( const int nr_in, const int nc_in, const bool flag_zero ) -{ - nr = nr_in; nc = nc_in; - c = new double[nr*nc]; - if(flag_zero) - memset( c, 0, sizeof(double)*nr*nc ); -} - -inline matrix Matrix_Wrapper::to_matrix() -{ - matrix m; - m.nr = nr; m.nc = nc; - m.c = c; - nr = nc = 0; - c = nullptr; - return m; -} - -} -#endif diff --git a/source/source_base/mcd.c b/source/source_base/mcd.c deleted file mode 100644 index 3ec039270f5..00000000000 --- a/source/source_base/mcd.c +++ /dev/null @@ -1,884 +0,0 @@ -#define MCD_VERBOSE -#define MCD_FASTFREE - -/* - * Softpixel MemCheckDeluxe v1.2.2 - * 2003 cwright - */ - -/* Copyright (C) 2002, 2003 Softpixel (http://softpixel.com/) - * This file is covered by a BSD-Style License, available in - * the LICENSE file in the root directory of this package, as well as at - * http://prj.softpixel.com/licenses/#bsd - */ - -#ifndef _MCD_GONE //if mcd is gone, don't do anything - -#ifdef _MCD_CHECK -#undef _MCD_CHECK //don't want trace allocations in this file =) -#endif - -#include -#ifndef __USE_GNU -#define __USE_GNU // activate strndup in GNU headers -#endif - -#include -#include -#include -//#include // fixes the vscanf warnings, - // but kills the va_* macros... - -#define MCD_C_FILE_BUILD //flag some stuff in mcd.h -#include "mcd.h" - -#ifdef WIN32 - #include -#else - #include - #include -#endif - -#define CHUNK_MAGIC 0x1337C0DE - -#ifndef MCD_FASTFREE -#define TOP_TO_BOTTOM 0x00 //search like old style -#define BOTTOM_TO_TOP 0xff //search from LastChunk up - -#define TRACK_LENGTH 5 //how many previous biases to use -#endif - -/* Number of bytes to add to each allocation to check for padding */ -//#define OVERFLOW_PAD 4 // unimplimented - -/* all allocations are accounted for with 'chunnks' */ -typedef struct ChunkS -{ -#ifdef MCD_FASTFREE - int magic; //help make sure the chunk is really a chunk -#endif - int size; //size of allocation in bytes - int line; //line it was alocated on -#ifdef WIN32 - LONGLONG id; // 64bit allocation ID -#else - long long id; // 64bit allocation ID -#endif - char * function = nullptr; //creating function - char * file = nullptr; //file function is in - void * ptr = nullptr; //pointer to allocation - struct ChunkS *next, //next chunk (null if nonw) - *prev; //previous chunk (null if nonw) -}Chunk; -/*note: *ptr+SOC is actually the allocation. with fastfree, ptr points to -the chunk itself, and SOC adds the offset where the actual data is */ - - -/* MemoryChunks stores the first entry to the allocation -chunks. Lastchunk stores the last one, to save list traversals. */ - -Chunk *MemoryChunks=NULL,*LastChunk=NULL; - -/* - smallest = smallest allocation - largest = largest allocation - current = amount currently allocated (for peak tracking) -*/ -int smallest=0,largest=0,maximum=0,current=0; - -/* - totalAlloc = total allocations - totalFree = total frees - peakAlloc = maximum number of outstanging allocations - totalOverflow = number of pointers that were overflowed - (unimplimented) -*/ -int totalAlloc=0, totalFree=0, peakAlloc=0, totalOverflow=0; - -/* Our log files (Null causes default behaviour) */ -FILE *RealTimeLog=NULL; -FILE *MemStatLog=NULL; - -//return realtimelog if it's set, otherwise stderr -#define RTL (RealTimeLog?RealTimeLog:stderr) - -/* SOC (SizeOf Chunk) is used in verbose mode to print correct pointers */ -#ifdef MCD_FASTFREE -#define SOC sizeof(Chunk) //fastfree pointer offset -#else -#define SOC 0 //no fastfree implies raw pointers (no offset) -#endif - -/* idbase is the id variable. each new allocation increments it. - to keep it useful, it is defined to be as long as possible - (64 bits) to keep it from wrapping (unless your program does - a _lot_ of allocating :) - */ -#ifdef WIN32 -LONGLONG idbase=0; -#else -long long idbase=0; -#endif - -/* chunklock is our simple lock variable to keep the chunk list thread-safe */ -char chunkLock=0; - -/* TrackBias is used to record the weight of the previous locations - TrackPos keeps our current position. each free increments - and wraps TrackPos. -*/ -#ifndef MCD_FASTFREE -unsigned int TrackBias[TRACK_LENGTH]={0x7f}; //start in middle -int TrackPos=0; -#endif - -/* Here, *p is the raw pointer. we stick our chunk pointer there for easy lookups */ -void addChunk(void*p,int size,char*fun,char*file,int line) -{ -#ifdef MCD_FASTFREE - Chunk *c=(Chunk*)p; //extra space for chunk is accounted for -#else //not FASTFREE - Chunk *c=(Chunk*)malloc(sizeof(Chunk)); -#endif - while(chunkLock) -#ifdef WIN32 - Sleep(1); -#else - usleep(1); -#endif - ++chunkLock; - #ifdef MCD_FASTFREE - c->magic=CHUNK_MAGIC; - #endif - c->size=size; - c->line=line; - c->function=fun; - c->file=file; - c->ptr=p; - c->next=NULL; - c->prev=NULL; - c->id=idbase++; - -//overflow padding checks (incomplete) -// pad=p+size; -// for(po=0;pomaximum) - maximum=current; - if(size>largest) - largest=size; - if(sizepeakAlloc) - peakAlloc=totalAlloc-totalFree; - - if(!MemoryChunks) //first chunk - { - MemoryChunks=c; - LastChunk=c; - } - else //new chunk - { - c->prev=LastChunk; - LastChunk->next=c; - LastChunk=c; - } - --chunkLock; - return; -} - -/* here, *p is the pointer, adjusted to be the chunk*/ -int delChunk(void*p/*,char*fun,char*file,int line*/) -{ - Chunk *s; -#ifndef MCD_FASTFREE - unsigned int bias,i; - int depth=0; -#endif - //should be handled in MCD_free now - //if(!p) - // return 0; // null pointer. list walker would crash - - while(chunkLock) -#ifdef WIN32 - Sleep(1); -#else - usleep(1); -#endif - ++chunkLock; //0wn the list - -#ifdef MCD_FASTFREE - - s=(Chunk*)p; - if(s->magic!=CHUNK_MAGIC) - { - /*freeing illegal pointer, or wrote out of bounds*/ - fprintf(RTL,"** Magic is [%x], which is wrong...[ptr:%x]\n",s->magic,(int)p); - fprintf(RTL,"** Either some buffer overflowed, or this is a bad pointer\n"); - fprintf(RTL,"** I'm just going to free it and see what happens... might segfault.\n"); - fflush(RTL); - /*because this is broken, we're not really sure what to do - with it... so, we opt to just free it*/ - --chunkLock; - return 1; - } - //pointer is to previous chunk, so use that as a basis, plus some pointer acrobatics - - if(s->prev==NULL) //first chunk - { - if(MemoryChunks->next) //next chunk gets promoted to MemoryChunks - { - current-=MemoryChunks->size; - //prev=MemoryChunks; - MemoryChunks=MemoryChunks->next; - MemoryChunks->prev=NULL; - } - else //only chunk left (empty list) - { - current-=MemoryChunks->size; - MemoryChunks=NULL; - LastChunk=NULL; - } - } - else if(s->next==NULL) //not first, is it last? - { - current-=s->size; - LastChunk=s->prev; - LastChunk->next=NULL; //last - } - else //must be in the middle somewhere... - { - - s->prev->next=s->next; - s->next->prev=s->prev; - current-=s->size; - } - ++totalFree; -#else //not FASTFREE - //non fastfree version. here, we expect the chunk to be freed upon exit. - - for(bias=0,i=0;iptr==p) - { - current-=s->size; - - if(s==MemoryChunks) //first in list - { - MemoryChunks=MemoryChunks->next; //set next as first - free(s); - if(MemoryChunks) //if there is a next, set prev to none - MemoryChunks->prev=NULL; - --chunkLock; - ++totalFree; - return 0; - } - else if(s==LastChunk) //last guy... - { - LastChunk=LastChunk->prev; - free(s); - if(LastChunk) //LastChunk could be first as well (no next), dont assume anything. - LastChunk->next=NULL; - --chunkLock; - ++totalFree; - return 0; - } - else // somewhere in the middle... - { - // s should have a next and first (handled by previous cases) - s->prev->next=s->next; - s->next->prev=s->prev; - free(s); - --chunkLock; - ++totalFree; - return 0; - } - } - else - { - if(bias==BOTTOM_TO_TOP) - s=s->prev; - else - s=s->next; - } - } - //hmm, didnt find it in the list, must be foreign or error... - #ifdef MCD_VERBOSE - --chunkLock; - fprintf(RTL,"** Pointer [%08X] wasn't in list...\n",(int)p); - fprintf(RTL,"** Either it came from an external library, or it's a bad pointer\n"); - fprintf(RTL,"** I'm just going to free it and see. might segfault\n"); - fflush(RTL); //mohan - #endif - //return 1; //error goes unnoticed for now -#endif //fastfree - --chunkLock; - return 0; -} - -/* (FastFree) because we use fake pointers, we must add sizeof(Chunk) to all of - them for proper pointer hiding. -*/ - -#ifndef WIN32 -inline //inline to make it a bit snappier -#endif -void phex(FILE*f, //print our 64bit id in a generic way -#ifdef WIN32 - LONGLONG v -#else - long long v -#endif -) -{ - int i; - char t[16]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'}; - for(i=0;i<64;i+=4) - { - fprintf(f,"%c",t[(v>>(60-i))&0xf]); - } - fprintf(f,"]\n"); - fflush(f); - return; -} - -void *MCD_malloc(int size,char*fun,char*file,int line) -{ -#ifdef MCD_FASTFREE - void*p=malloc(size+sizeof(Chunk)); -#else //not FASTFREE - void*p=malloc(size); -#endif - -#ifdef MCD_VERBOSE - fprintf(RTL,"malloc\t[%X][%i bytes]\tfrom %s:%s, line %i [id ", - (int)p+SOC,size,file,fun,line); - phex(RTL,idbase); -#endif - addChunk(p,size,fun,file,line); - - return p -#ifdef MCD_FASTFREE - +sizeof(Chunk) -#endif - ; -} - - -void *MCD_calloc(int size,char*fun,char*file,int line) -{ -#ifdef MCD_FASTFREE - void*p=calloc(size+sizeof(Chunk),1); -#else - void*p=calloc(size,1); -#endif -#ifdef MCD_VERBOSE - fprintf(RTL,"calloc\t[%X][%i bytes]\tfrom %s:%s, line %i [id ", - (int)p+SOC,size,file,fun,line); - phex(RTL,idbase); -#endif - - addChunk(p,size,fun,file,line); - return p -#ifdef MCD_FASTFREE - +sizeof(Chunk) -#endif - ; -} - -/* instead of really reallocing, we just allocate a new area, and copy - this shouldn't matter, but there may be a better way */ -void *MCD_realloc(void *p,int size,char*fun,char*file,int line) -{ -#ifdef MCD_FASTFREE - Chunk *New,*Old; - - void*n=malloc(size+sizeof(Chunk)); - New=(Chunk*)n; - Old=(Chunk*)p-sizeof(Chunk); - - if(Old->size==size) //no change?? - return p; - - if(size>Old->size) //allocated more room, just copy old stuff - memcpy(n,p-sizeof(Chunk),Old->size+sizeof(Chunk)); - else //allocated less room (truncate) - memcpy(n,p-sizeof(Chunk),Old->size+sizeof(Chunk)); - - -#ifdef MCD_VERBOSE - fprintf(RTL,"realloc\t[%X]->[%X][%i bytes]\tfrom %s:%s, line %i [id ", - (int)p+SOC,(int)n+SOC,size,file,fun,line); - phex(RTL,idbase); -#endif - - delChunk(p/*,fun,file,line*/); - free(p-sizeof(Chunk)); - addChunk(n,size,fun,file,line); - return n+sizeof(Chunk); -#else //non-fastfree - void*n=realloc(p,size); -#ifdef MCD_VERBOSE - fprintf(RTL,"realloc [%X]->[%X][%ib] from %s:%s, line %i [id ", - (int)p+SOC,(int)n+SOC,size,file,fun,line); - phex(RTL,idbase); -#endif - delChunk(p/*,fun,file,line*/); - addChunk(n,size,fun,file,line); - return n; -#endif -} - -// getcwd isn't implimented in fastfree yet -#ifndef MCD_FASTFREE -#ifndef WIN32 //win32's getcwd doesn't allocate anyway -char *MCD_getcwd(char *p,int size, char*fun,char*file,int line) -{ - int olderrno,wsize=256; - char *tmp=NULL; - - //according to man page: - // if p is null, size bytes are allocated - // if p is null and size is zero, strlen(wd) bytes are allocated - - if(p!=NULL) //ok, do nothing - { - return getcwd(p,size); - } - else // ok, interesting stuff - { - if(!size) //alloc as big as needed - { - while(tmp==NULL) - { - tmp=(char*)malloc(wsize); - olderrno=errno; - if(getcwd(tmp,0)==NULL) - { //too small, let's try again - wsize*=2; - free(tmp); - } - errno=olderrno; //restore old errno, just in case - } - p=(char*)malloc(strlen(tmp)+1); - memcpy(p,tmp,strlen(tmp)+1); - free(tmp); - addChunk(p, strlen(tmp)+1, file, fun, line); - } - else //alloc size bytes - { - // TODO: see if failure condition here still allocs p - p=(char*)malloc(size); - addChunk(p, size, file, fun, line); - getcwd(p,size); - } - return p; - } -} -#endif //win32 -#endif //fastfree - -/*(in fastfree) can't make this simple due to pointer stuffing :( */ - -#ifdef __GNUC__ -char *MCD_strdup(const char*s,char*fun,char*file,int line) -#else -char *MCD_strdup(char*s,char*fun,char*file,int line) -#endif -{ - char * n = nullptr; - -#ifdef MCD_FASTFREE - n=(char*)malloc(sizeof(char)*strlen(s)+1+sizeof(Chunk)); - memcpy(n+sizeof(Chunk),s,strlen(s)); -#else //non-fastfree - n=strdup(s); //unpointer'd version -#endif -#ifdef MCD_VERBOSE - fprintf(RTL,"strdup\t[%X][%ib]\t\tfrom %s:%s, line %i [id ", - (int)n+SOC,(int)strlen(s),file,fun,line); - phex(RTL,idbase); -#endif - addChunk(n,strlen(s)+1,fun,file,line); - return n -#ifdef MCD_FASTFREE - +sizeof(Chunk) -#endif - ; -} - -#ifndef WIN32 //win32 lacks strndup completely -#ifdef __GNUC__ -char *MCD_strndup(const char*s, int z,char*fun,char*file,int line) -#else -char *MCD_strndup(char*s, int z,char*fun,char*file,int line) -#endif -{ - char * n = nullptr; - int size; - if((signed)strlen(s)>z) - size=z; - else - size=strlen(s); -#ifdef MCD_FASTFREE - n=(char*)malloc(sizeof(char)*size+1+sizeof(Chunk)); - memcpy(n+sizeof(Chunk),s,size); -#else //non-fastfree - n=(char*)strndup(s,z); //unpointer'd version -#endif -#ifdef MCD_VERBOSE - fprintf(RTL,"strndup\t[%X][%ib]\t\tfrom %s:%s, line %i [id ", - (int)n+SOC,(int)size,file,fun,line); - phex(RTL,idbase); -#endif - addChunk(n,size+1,fun,file,line); - return n -#ifdef MCD_FASTFREE - +sizeof(Chunk) -#endif - ; -} - -/* asprintf/vasprintf based on a patch from Stephen Lee */ -#ifdef _GNU_SOURCE - -int MCD_asprintf(char **ptr,const char *fmt,char *fun, char*file, int line,...) -{ - int retval; - va_list argptr; -#ifdef MCD_FASTFREE - void * fc = nullptr; -#endif - - va_start(argptr,line); - if((retval=vasprintf(ptr, fmt, argptr))<0) - { - #ifdef MCD_VERBOSE - fprintf(RTL,"asprintf failure %s:%s, line %i [id ", - file,fun,line); - #endif - return retval; - } - va_end(argptr); - -#ifdef MCD_FASTFREE - fc=malloc(sizeof(Ghunk)+retval+1); //make string+chunk - memcpy(fc+sizeof(Chunk),ptr,retval); //hack in string - free(ptr); //clean up string - ptr=fc+sizeof(Chunk); //hack pointer - addChunk(fc,retval+1,fun,file,line); -#else //non-fastfree - addChunk(*ptr, retval+1, fun, file, line); -#endif - #ifdef MCD_VERBOSE - fprintf(RTL,"asprintf\t[%X][%ib]\t\tfrom %s:%s, line %i [id ", - (int)*ptr+SOC,retval+1,file,fun,line); - phex(RTL,idbase); - #endif - return retval; -} - -int MCD_vasprintf(char **ptr,const char *fmt,va_list argptr,char *fun, char*file, int line) -{ - int retval; -#ifdef MCD_FASTFREE - void * fc = nullptr; -#endif - - if((retval=vasprintf(ptr, fmt, argptr))<0) - { - #ifdef MCD_VERBOSE - fprintf(RTL,"vasprintf failure %s:%s, line %i [id ", - file,fun,line); - #endif - return retval; - } -#ifdef MCD_FASTFREE - fc=malloc(sizeof(Chunk)+retval+1); - memcpy(fc+sizeof(Chunk),ptr,retval); - free(ptr); - ptr=fc+sizeof(Chunk); - addChunk(fc,retval+1,fun,file,line); -#else //non-fastfree - addChunk(*ptr, retval+1, fun, file, line); -#endif - #ifdef MCD_VERBOSE - fprintf(RTL,"vasprintf\t[%X][%ib]\t\tfrom %s:%s, line %i [id ", - (int)*ptr+SOC,retval+1,file,fun,line); - phex(RTL,idbase); - #endif - return retval; -} - -#endif // _GNU_SOURCE -#endif // WIN32 - -/* scanf family patch provided by Stephen Lee */ -void scan_args(const char*fmt,va_list argptr,char*fun,char*file,int line); - -#ifndef __cplusplus // c++ doesn't seem to compile when vscanf etc are implicit -#ifndef WIN32 //win32 doesn't link vscanf stuff nicely -int MCD_scanf(const char *fmt,char*fun,char*file,int line,...) -{ - int retval; - va_list argptr; - - va_start(argptr,line); - if((retval=vscanf(fmt,argptr))<1) { - /* no args, so let's return */ - va_end(argptr); - return retval; - } - va_end(argptr); - va_start(argptr,line); - scan_args(fmt,argptr,fun,file,line); - va_end(argptr); - return retval; -} - -int MCD_fscanf(FILE *stream,const char *fmt,char*fun,char*file,int line,...) -{ - int retval; - va_list argptr; - - va_start(argptr,line); - if((retval=vfscanf(stream,fmt,argptr))<1) { - /* no args, so let's return */ - va_end(argptr); - return retval; - } - va_end(argptr); - va_start(argptr,line); - scan_args(fmt,argptr,fun,file,line); - va_end(argptr); - return retval; -} - -int MCD_sscanf(const char *str,const char *fmt,char*fun,char*file,int line,...) -{ - int retval; - va_list argptr; - - va_start(argptr,line); - if((retval=vsscanf(str,fmt,argptr))<1) { - /* no args, so let's return */ - va_end(argptr); - return retval; - } - va_end(argptr); - va_start(argptr,line); - scan_args(fmt,argptr,fun,file,line); - va_end(argptr); - return retval; -} -#endif //Win32 link -#endif //__cplusplus (implicit blah blah) - -/* scanf etc helper function */ -void scan_args(const char *fmt,va_list argptr,char*fun,char*file,int line) -{ - char ** ptr = nullptr; - void * dummy = nullptr; // clear up the unused warning - - for(;*fmt;fmt++) { - if(*fmt!='%') - continue; - switch(*(++fmt)) { - case 'a': /* malloc'd string */ - ptr=(char **)va_arg(argptr,char *); - addChunk(*ptr,strlen(*ptr)+1,fun,file,line); - break; - case '%': - break; - default: /* next arg */ - dummy=va_arg(argptr,void *); - break; - } - } - return; -} - -void MCD_free(void *p,char*fun,char*file,int line) -{ - //quiet the unused argument warnings when MCD_VERBOSE isn't defined -#ifndef MCD_VERBOSE - void *dummy=&fun; - dummy=&file; - dummy=(char*)&line; -#endif //MCD_VERBOSE - - if(!p) - { -#ifdef MCD_VERBOSE - fprintf(RTL,"** Call to free with NULL argument in %s:%s, line %i\n", - file,fun,line); -#endif //mcd_verbose - /* If they want native free(null), do it */ - /* in case free(null) does something bad, */ - /* rather than let them think it's ok */ - /* Must be enabled, off by default. */ - #ifdef MCD_FREE_NULL - free(p); //native free(null) for same behaviour - #endif - return; //we do nothing - } -#ifdef MCD_VERBOSE - fprintf(RTL,"free\t[%X]\t\tfrom %s:%s, line %i\n", - (int)p,file,fun,line); -#endif - //should find id, but that's delchunk's job... - -#ifdef MCD_FASTFREE - // if delChunk returns something, it means the pointer is crazy - if(!delChunk(p-sizeof(Chunk)/*,file,fun,line*/)) - free(p-sizeof(Chunk)); - else - { - fprintf(RTL,"** Bad or Foreign pointer [%08X] from %s:%s line %i\n", - (int)p,file,fun,line); - fflush(RTL); - free(p); //assume it's just a foreign pointer (could die) - } -#else //not fastfree - - delChunk(p); //return value here doesn't matter, but 1 is still 'error' - // error behaviour is simply the same (free and see) - free(p); -#endif - return; -} - -void showMemStats() -{ - Chunk*c=MemoryChunks; - int total=0; - - FILE * o = nullptr; - - if(MemStatLog) - o=MemStatLog; - else - o=stdout; - - fprintf(o,"Memory Stats:\n"); - while(c) - { - fprintf(o,"%12i bytes allocated in [%s:%s], line %7i [%08X][id ", - c->size,c->function,c->file,c->line,(int)c->ptr+SOC); - phex(o,c->id); - total+=c->size; - c=c->next; - fflush(o); - } - fprintf(o,"\t\tTotal unfreed bytes: %12i\n",total); - //fprintf(o,"\t\tOverflows[approx.]: %12i\n",totalOverflow); - fprintf(o,"\t\tPeak memory usage: %12i\n",maximum); - fprintf(o,"\t\tLargest allocation: %12i\n",largest); - fprintf(o,"\t\tSmallest allocation: %12i\n",smallest); - fprintf(o,"\t\tTotal Allocations: %12i\n",totalAlloc); - fprintf(o,"\t\tTotal Frees: %12i\n",totalFree); - fprintf(o,"\t\tPeak Allocations: %12i\n",peakAlloc); - fprintf(o,"\n"); - return; -} - -void _MCD_RealTimeLog(FILE*f) -{ - RealTimeLog=f; - return; -} - -void _MCD_MemStatLog(FILE*f) -{ - MemStatLog=f; - return; -} - -#ifdef __cplusplus - -char *_MCD_LastSetFile,*_MCD_LastSetFun; -int _MCD_LastSetLine; - - -/*void* operator new (unsigned int size,char *file, char*fun) -{ - return MCD_malloc(size,file,fun,5); -} - -void* operator new[] (unsigned int size,char*file, char*fun) -{ - return MCD_malloc(size,file,fun,5); -}*/ - - -/*void* operator new (size_t size, char* file, char*fun, int line) -{ - return MCD_malloc(size, file, fun, line); -} - -void* operator new[] (size_t size, char* file, char* fun, int line) -{ - return MCD_malloc(size, file, fun, line); -} - -void operator delete(void* buffer, char*file, char*fun, int line) -{ - MCD_free(buffer, file, fun, line); - -void operator delete[] (void* buffer, char*file, char*fun, int line) -{ - MCD_free(buffer, file, fun, line); - -}*/ - -/* -void* operator new (size_t size) -{ - return MCD_malloc(size, __FILE__, __FUNCTION__, __LINE__); -} - -void* operator new[] (size_t size) -{ - return MCD_malloc(size, __FILE__, __FUNCTION__, __LINE__); -} - -void operator delete(void* buffer) -{ - MCD_free(buffer, __FILE__, __FUNCTION__, __LINE__); -} - -void operator delete[](void* buffer) -{ - MCD_free(buffer, __FILE__, __FUNCTION__, __LINE__); -}*/ -#endif // __cplusplus - -#endif //_MCD_GONE diff --git a/source/source_base/mcd.h b/source/source_base/mcd.h deleted file mode 100644 index 494c2f7e46f..00000000000 --- a/source/source_base/mcd.h +++ /dev/null @@ -1,285 +0,0 @@ -/* - * Softpixel MemCheckDeluxe v1.2.2 - * 2003 cwright - */ - -/* Copyright (C) 2002, 2003 Softpixel (http://softpixel.com/) - * This file is covered by a BSD-Style License, available in - * the LICENSE file in the root directory of this package, as well as at - * http://prj.softpixel.com/licenses/#bsd - */ - - -#ifndef MCD_H -#define MCD_H - -#include //for file redirection - -#define MCD_VERSION 0x010201 - -#ifdef __cplusplus -extern "C" { -#endif - -#ifndef _MCD_GONE - -#ifdef _MCD_CHECK - -#include - -#ifdef WIN32 - #define __FUNCTION__ __FILE__ -#endif - -#ifndef MCD_QUIET - #ifndef WIN32 - #warning - Building with memory checking. - #warning expect lower performance. - - #endif -#endif - -//warn about redefinitions (sometimes these can be #defined) - -#ifdef malloc -#undef malloc - #ifndef WIN32 - #warning ------ Redefining malloc() ------ - #endif -#endif - -#ifdef calloc -#undef calloc - #ifndef WIN32 - #warning ------ Redefining calloc() ------ - #endif -#endif - -#ifdef realloc -#undef realloc - #ifndef WIN32 - #warning ------ Redefining realloc() ------ - #endif -#endif - -#ifdef free -#undef free - #ifndef WIN32 - #warning ------ Redefining free() ------ - #endif -#endif - -#ifdef strdup -#undef strdup - #ifndef WIN32 - #warning ------- Redefining strdup() -------- - #endif -#endif - -#ifdef strndup -#undef strndup - #ifndef WIN32 - #warning ------- Redefining strndup() ------- - #endif -#endif - -#ifdef asprintf -#undef asprintf - #ifndef WIN32 - #warning ------ Redefining asprintf() ------ - #endif -#endif - -#ifdef vasprintf -#undef vasprintf - #ifndef WIN32 - #warning ------ Redefining vasprintf() ------ - #endif -#endif - -#ifdef scanf -#undef scanf - #ifndef WIN32 - #warning ------ Redefining scanf() ------ - #endif -#endif - -#ifdef fscanf -#undef fscanf - #ifndef WIN32 - #warning ------ Redefining fscanf() ------ - #endif -#endif - -#ifdef sscanf -#undef sscanf - #ifndef WIN32 - #warning ------ Redefining sscanf() ------ - #endif -#endif - -#ifdef getcwd -#undef getcwd - #ifndef WIN32 - #warning ------ Redefining getcwd() ------ - #endif -#endif - -#define strdup(p) MCD_strdup(p,__FUNCTION__,__FILE__,__LINE__) -#ifndef WIN32 //no strndup in win32 -#define strndup(p,n) MCD_strndup(p,n,__FUNCTION__,__FILE__,__LINE__) -#endif -#define malloc(size) MCD_malloc(size,__FUNCTION__,__FILE__,__LINE__) -#define calloc(n,s) MCD_calloc(s*n,__FUNCTION__,__FILE__,__LINE__) -#define realloc(p,s) MCD_realloc(p,s,__FUNCTION__,__FILE__,__LINE__) - -#ifdef _GNU_SOURCE -#define asprintf(p,f,args...) MCD_asprintf(p,f,__FUNCTION__,__FILE__,__LINE__, ## args) -#define vasprintf(p,f,ap) MCD_vasprintf(p,f,ap,__FUNCTION__,__FILE__,__LINE__) -#endif - -#ifndef WIN32 -/* windows doesn't like variable arguments in #define's */ -#ifndef MCD_FASTFREE //fastfree isnt supported here yet -#define scanf(f,args...) MCD_scanf(f,__FUNCTION__,__FILE__,__LINE__, ## args) -#define fscanf(s,f,args...) MCD_fscanf(s,f,__FUNCTION__,__FILE__,__LINE__, ## args) -#define sscanf(s,f,args...) MCD_sscanf(s,f,__FUNCTION__,__FILE__,__LINE__, ## args) -#define getcwd(p,s) MCD_getcwd(p,s,__FUNCTION__,__FILE__,__LINE__) -#endif //fastfree -#endif //win32 - -#define free(p) MCD_free(p,__FUNCTION__,__FILE__,__LINE__) - -#endif // _MCD_CHECK - - -/* !!! These are called by the defines only. Do NOT use directly. !!! */ -void *MCD_malloc(int size, char*, char*, int); -void *MCD_calloc(int size, char*, char*, int); -void *MCD_realloc(void *p, int size, char*, char*, int); -char *MCD_getcwd (char *p, int size, char*, char*, int); -#ifdef __GNUC__ -char *MCD_strdup (const char*s, char*, char*, int); -char *MCD_strndup(const char*s, int n, char*, char*, int); -#else -char *MCD_strdup (char*s, char*, char*, int); -char *MCD_strndup(char*s, int n, char*, char*, int); -#endif -int MCD_scanf (const char *fmt, char*fun, char*file, int line,...); -int MCD_fscanf(FILE *stream,const char *fmt,char*fun,char*file,int line,...); -int MCD_sscanf(const char *str,const char *fmt,char*fun,char*file,int line,...); -// Private MCD function, no need to be in public namespace -//void scan_args(const char*fmt,va_list argptr,char*fun,char*file,int line); -#ifdef _GNU_SOURCE -int MCD_asprintf(char **ptr, const char*fmt,char*,char*,int,...); -int MCD_vasprintf(char **ptr, const char*fmt,va_list argptr,char*,char*,int); -#endif // _GNU_SOURCE - -void MCD_free(void *p,char*,char*,int); - -/* --- call this for memory stats --- */ -void showMemStats(void); - -/* --- to send realtime stats somewhere other than stderr, - put an opened fp in here --- */ -void _MCD_RealTimeLog(FILE*); - -/* --- to send showMemStats() somewhere other than stdout, - put an opened fp in here --- */ -void _MCD_MemStatLog(FILE*); - -#else //MCD is gone - -/* define functions so source will compile without modification */ -#define showMemStats() -#define _MCD_RealTimeLog(x) -#define _MCD_MemStatLog(x) - -#endif //_MCD_GONE - -#ifdef __cplusplus -} /* extern "C" */ -#endif - -#ifndef _MCD_GONE // check again, outside of extern "C" mode -#ifdef __cplusplus /* Some C++ new/delete operator overloading */ - -#ifndef WIN32 - #warning C++ Extentions Enabled -#endif - -#ifdef new - #undef new -#endif - -#ifdef delete - #undef delete -#endif - -extern char *_MCD_LastSetFile,*_MCD_LastSetFun; -extern int _MCD_LastSetLine; - -inline void setFileFunLineState(char*file,char*fun,int line) -{ - _MCD_LastSetLine=line; - _MCD_LastSetFile=file; - _MCD_LastSetFun=fun; -} - -inline void* operator new (unsigned int size,char *file, -char*fun,int line) -{ - return MCD_malloc(size,file,fun,line); -} - -inline void* operator new[] (unsigned int size,char*file, -char*fun,int line) -{ - return MCD_malloc(size,file,fun,line); -} - -// currently, passing args to delete operator is not working at all... -#ifndef WIN32 //win32 doesn't like default params to delete -inline void operator delete (void * buf,char*file=__FILE__, -#ifdef WIN32 -char*fun=__FILE__,int line=__LINE__) -#else -char*fun=__FUNCTION__,int line=__LINE__) -#endif -{ - MCD_free(buf,_MCD_LastSetFile,_MCD_LastSetFun,_MCD_LastSetLine); -} -// ...so we have these here for the day they work, which isn't today. -inline void operator delete[] (void * buf,char*file=__FILE__, -#ifdef WIN32 -char*fun=__FILE__,int line=__LINE__) -#else -char*fun=__FUNCTION__,int line=__LINE__) -#endif -{ - MCD_free(buf,_MCD_LastSetFile,_MCD_LastSetFun,_MCD_LastSetLine); -} -#endif // win32 default delete params - -inline void operator delete (void * buf) -{ - MCD_free(buf,0,0,0); -} - -inline void operator delete[] (void * buf) -{ - MCD_free(buf,0,0,0); -} - -#ifdef WIN32 -#define new new(__FILE__,__FILE__,__LINE__) -#else -#define new new(__FILE__,__FUNCTION__,__LINE__) -#endif -//#define delete setFileFunLineState(__FILE__,__FUNCTION__,__LINE__);delete -//#define delete delete(__FILE__,__FUNCTION__,__LINE__) - -#endif // __cplusplus -#endif // _MCD_GONE - - -#endif // MCD_H diff --git a/source/source_base/test/CMakeLists.txt b/source/source_base/test/CMakeLists.txt index ce9df8138ba..16ba949d640 100644 --- a/source/source_base/test/CMakeLists.txt +++ b/source/source_base/test/CMakeLists.txt @@ -218,12 +218,6 @@ AddTest( LIBS parameter base device ) -AddTest( - TARGET MODULE_BASE_projgen - SOURCES projgen_test.cpp - LIBS parameter base device -) - AddTest( TARGET MODULE_BASE_clebsch_gordan_coeff_test SOURCES cg_coeff_test.cpp diff --git a/source/source_base/vector3.h b/source/source_base/vector3.h index 5fb2bf3f7e1..27a7b60d412 100644 --- a/source/source_base/vector3.h +++ b/source/source_base/vector3.h @@ -7,9 +7,6 @@ #include #include -#ifdef _MCD_CHECK -#endif - namespace ModuleBase { // Small epsilon value for numerical comparisons diff --git a/source/source_basis/module_nao/CMakeLists.txt b/source/source_basis/module_nao/CMakeLists.txt index 29e091510f3..a5fd92bb193 100644 --- a/source/source_basis/module_nao/CMakeLists.txt +++ b/source/source_basis/module_nao/CMakeLists.txt @@ -1,3 +1,16 @@ +# smoothgen is consumed by both NAO (atomic_radials) and PW +# (source_pw/module_proj/onsite_proj_init) code, so it must be built +# and linked regardless of ENABLE_LCAO. +add_library( + nao_smoothgen + OBJECT + smoothgen.cpp +) + +if(ENABLE_COVERAGE) + add_coverage(nao_smoothgen) +endif() + if(ENABLE_LCAO) add_library( numerical_atomic_orbitals @@ -21,6 +34,7 @@ if(ENABLE_LCAO) endif() if(BUILD_TESTING) + add_subdirectory(unittests) if(ENABLE_MPI) add_subdirectory(test) endif() diff --git a/source/source_basis/module_nao/atomic_radials.cpp b/source/source_basis/module_nao/atomic_radials.cpp index edcf1161a43..3827003893f 100644 --- a/source/source_basis/module_nao/atomic_radials.cpp +++ b/source/source_basis/module_nao/atomic_radials.cpp @@ -7,7 +7,7 @@ // FIXME: should update with pyabacus // #include "source_io/orb_io.h" -#include "source_base/projgen.h" +#include "source_basis/module_nao/smoothgen.h" #include #include @@ -102,7 +102,6 @@ void AtomicRadials::build(RadialSet* const other, const int itype, const double std::vector rvalue_new; smoothgen(ngrid, rgrid, rvalue, rcut, rvalue_new); ngrid = rvalue_new.size(); - // projgen(l, ngrid, rgrid, rvalue, rcut, 20, rvalue_new); // build the new on-site orbitals this->chi_[ichi].build(l, true, ngrid, rgrid, rvalue_new.data(), 0, izeta, symbol_, itype, false); } diff --git a/source/source_base/projgen.cpp b/source/source_basis/module_nao/smoothgen.cpp similarity index 65% rename from source/source_base/projgen.cpp rename to source/source_basis/module_nao/smoothgen.cpp index 681bc2a6f06..42d9dd50e8a 100644 --- a/source/source_base/projgen.cpp +++ b/source/source_basis/module_nao/smoothgen.cpp @@ -1,84 +1,16 @@ -#include "projgen.h" +#include "smoothgen.h" #include #include #include -#include #include +#include #include "source_base/cubic_spline.h" #include "source_base/math_integral.h" -#include "source_base/math_sphbes.h" - -using namespace ModuleBase; - -void projgen(const int l, - const int nr, - const double* r, - const double* chi, - const double rcut, - const int nbes, - std::vector& alpha) -{ - assert(rcut < r[nr - 1]); - assert(std::is_sorted(r, r + nr)); - - std::vector dr(nr); - std::adjacent_difference(r, r + nr, dr.begin()); - - // lower_bound returns the first element that is equal or larger than rcut - int nr_proj = std::distance(r, std::lower_bound(r, r + nr, rcut)) + 1; - - // zeros of spherical Bessel function - std::vector theta(nbes); - Sphbes::sphbes_zeros(l, nbes, theta.data()); - - // z & w vectors (see notes) - std::vector z(nbes); - std::vector w(nbes); - - std::transform(theta.begin(), theta.end(), z.begin(), [rcut, l](double theta_p) { - return 0.5 * std::pow(rcut, 3) * std::pow(Sphbes::sphbesj(l + 1, theta_p), 2); - }); - - // r^2 * chi (independent from p) - std::vector tmp(nr_proj); - std::transform(r, r + nr_proj, chi, tmp.begin(), [](double r_i, double chi_i) { return r_i * r_i * chi_i; }); - - // r^2 * chi * j_l(theta[p] * r / rcut) (dependent on p) - std::vector integrand(nr_proj); - - for (int p = 0; p < nbes; ++p) - { - std::transform(r, r + nr_proj, tmp.begin(), integrand.begin(), [theta, p, rcut, l](double r_i, double tmp_i) { - return tmp_i * Sphbes::sphbesj(l, theta[p] * r_i / rcut); - }); - w[p] = Integral::simpson(nr_proj, integrand.data(), &dr[1]); - } - - // optimal coefficients - std::vector c(nbes, 0.0); - std::transform(w.begin(), w.end(), z.begin(), c.begin(), [](double w_p, double z_p) { return w_p * w_p / z_p; }); - double prefac = 1.0 / std::sqrt(std::accumulate(c.begin(), c.end(), 0.0)); - std::transform(w.begin(), w.end(), z.begin(), c.begin(), [prefac](double w_p, double z_p) { - return prefac * w_p / z_p; - }); - - // new radial function - alpha.resize(nr_proj); - std::fill(alpha.begin(), alpha.end(), 0.0); - for (int i = 0; i < nr_proj; ++i) - { - for (int p = 0; p < nbes; ++p) - { - alpha[i] += c[p] * Sphbes::sphbesj(l, theta[p] * r[i] / rcut); - } - } -} void smoothgen(const int nr, const double* r, const double* chi, const double rcut, std::vector& alpha) { - // lambda function for generate the new radial function assert(rcut < r[nr - 1]); assert(std::is_sorted(r, r + nr)); @@ -127,7 +59,9 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc std::vector dalpha(nr_proj); cubspl_alpha.eval(nr_proj, r, nullptr, dalpha.data()); for (int i = 0; i < nr_proj; i++) + { dalpha[i] -= dchi[i]; + } // r^2 * dchi (independent from p) std::vector tmp(nr_proj); std::transform(r, r + nr_proj, dalpha.data(), tmp.begin(), [](double r_i, double dalpha_i) { @@ -141,7 +75,7 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc dalpha.data() + nr_proj, tmp.begin(), integrand.begin(), - [](double dalpha_i, double tmp_i) { return tmp_i * dalpha_i; }); + [](double dalpha_i, double tmp_i) { return dalpha_i * tmp_i; }); return ModuleBase::Integral::simpson(nr_proj, integrand.data(), &dr[1]); }; diff --git a/source/source_basis/module_nao/smoothgen.h b/source/source_basis/module_nao/smoothgen.h new file mode 100644 index 00000000000..4c1c1832c67 --- /dev/null +++ b/source/source_basis/module_nao/smoothgen.h @@ -0,0 +1,29 @@ +#ifndef SMOOTHGEN_H_ +#define SMOOTHGEN_H_ + +#include + +/** + * @brief Smoothly truncates a radial function at a cutoff radius. + * + * Starting with a radial function chi on grid r, this function generates a + * new radial function alpha on the truncated grid (points below rcut) such + * that alpha decays smoothly to zero at rcut. The smoothness parameter sigma + * is optimized by bisection so that the overlap between the derivatives of + * alpha and chi is minimized, and alpha is normalized so that the overlap + * integral with chi equals one. + * + * @param[in] nr number of grid points + * @param[in] r radial grid + * @param[in] chi radial function + * @param[in] rcut cutoff radius + * @param[out] alpha smoothed radial function on the truncated grid + */ +void smoothgen( + const int nr, + const double* r, + const double* chi, + const double rcut, + std::vector& alpha); + +#endif diff --git a/source/source_basis/module_nao/test/CMakeLists.txt b/source/source_basis/module_nao/test/CMakeLists.txt index 626023b7d24..f011f71faeb 100644 --- a/source/source_basis/module_nao/test/CMakeLists.txt +++ b/source/source_basis/module_nao/test/CMakeLists.txt @@ -12,6 +12,7 @@ AddTest( SOURCES atomic_radials_test.cpp ../atomic_radials.cpp + ../smoothgen.cpp ../radial_set.cpp ../numerical_radial.cpp ../../module_ao/orb_atomic_lm.cpp @@ -73,6 +74,7 @@ AddTest( radial_collection_test.cpp ../radial_collection.cpp ../atomic_radials.cpp + ../smoothgen.cpp ../beta_radials.cpp ../hydrogen_radials.cpp ../pswfc_radials.cpp @@ -91,6 +93,7 @@ AddTest( ../two_center_table.cpp ../radial_collection.cpp ../atomic_radials.cpp + ../smoothgen.cpp ../beta_radials.cpp ../hydrogen_radials.cpp ../pswfc_radials.cpp @@ -122,6 +125,7 @@ AddTest( ../real_gaunt_table.cpp ../radial_collection.cpp ../atomic_radials.cpp + ../smoothgen.cpp ../beta_radials.cpp ../hydrogen_radials.cpp ../pswfc_radials.cpp @@ -143,6 +147,7 @@ AddTest( ../real_gaunt_table.cpp ../radial_collection.cpp ../atomic_radials.cpp + ../smoothgen.cpp ../beta_radials.cpp ../hydrogen_radials.cpp ../pswfc_radials.cpp @@ -156,13 +161,13 @@ AddTest( AddTest( TARGET MODULE_NAO_generate_projection SOURCES - projgen_test.cpp ../two_center_bundle.cpp ../two_center_integrator.cpp ../two_center_table.cpp ../real_gaunt_table.cpp ../radial_collection.cpp ../atomic_radials.cpp + ../smoothgen.cpp ../beta_radials.cpp ../hydrogen_radials.cpp ../pswfc_radials.cpp diff --git a/source/source_basis/module_nao/test/projgen_test.cpp b/source/source_basis/module_nao/test/projgen_test.cpp deleted file mode 100644 index 620bd7786d6..00000000000 --- a/source/source_basis/module_nao/test/projgen_test.cpp +++ /dev/null @@ -1,136 +0,0 @@ -#include "source_base/projgen.h" -#include "gtest/gtest.h" - -#include "source_base/math_integral.h" -#include "source_base/math_sphbes.h" -#include "source_base/cubic_spline.h" - -#include -#include -#include -#include -#include - -/*********************************************************** - * Unit test of functions in projgen.cpp - * "projgen" : generate the projector coefficients - * "smoothgen" : smooth the projector coefficients and optimize the sigma - * - ***********************************************************/ - -TEST(projgen_test, projgen) -{ - // test orbital r^2 * exp(-r) - int l = 2; - double dr = 0.01; - double rcut_nao = 10; - int nr_nao = int(rcut_nao / dr) + 1; - std::vector r(nr_nao); - std::vector orb(nr_nao); - - for (int i = 0; i < nr_nao; ++i) { - r[i] = i * dr; - orb[i] = r[i] * r[i] * std::exp(-r[i]); - } - - // normalize the input orbital - std::vector integrand(nr_nao); - std::transform(r.begin(), r.end(), orb.begin(), integrand.begin(), - [](double r_i, double orb_i) { return std::pow(r_i * orb_i, 2); }); - double N = 1.0 / std::sqrt(ModuleBase::Integral::simpson(nr_nao, integrand.data(), dr)); - std::for_each(orb.begin(), orb.end(), [N](double& chi_i) { chi_i *= N; }); - - // projector information - double rcut_proj = 7.0; - int nbes = 7; - std::vector alpha; - - projgen(l, nr_nao, r.data(), orb.data(), rcut_proj, nbes, alpha); - - // compare with python script result - std::vector ref = { - 0.000000000000e+00, - 2.344902364599e-05, - 9.378381332712e-05, - 2.109675345121e-04, - 3.749388271050e-04, - 5.856118515995e-04, - 8.428763536364e-04, - 1.146597746904e-03, - 1.496617214310e-03, - 1.892751827321e-03, - 2.334794683381e-03, - 2.822515061259e-03, - 3.355658594204e-03, - 3.933947460740e-03, - 4.557080592928e-03, - 5.224733901903e-03, - 5.936560520491e-03, - 6.692191062668e-03, - 7.491233899644e-03, - 8.333275452302e-03, - }; - - for (int i = 0; i < 20; ++i) { - EXPECT_NEAR(alpha[i], ref[i], 1e-12); - } -} - -TEST(smoothgen_test, smoothgen) -{ - // test orbital r^2 * exp(-r) - int l = 2; - double dr = 0.01; - double rcut_nao = 10; - int nr_nao = int(rcut_nao / dr) + 1; - std::vector r(nr_nao); - std::vector orb(nr_nao); - - for (int i = 0; i < nr_nao; ++i) { - r[i] = i * dr; - orb[i] = r[i] * r[i] * std::exp(-r[i]); - } - - // normalize the input orbital - std::vector integrand(nr_nao); - std::transform(r.begin(), r.end(), orb.begin(), integrand.begin(), - [](double r_i, double orb_i) { return std::pow(r_i * orb_i, 2); }); - double N = 1.0 / std::sqrt(ModuleBase::Integral::simpson(nr_nao, integrand.data(), dr)); - std::for_each(orb.begin(), orb.end(), [N](double& chi_i) { chi_i *= N; }); - - // projector information - double rcut_proj = 7.0; - int nbes = 7; - std::vector alpha; - - smoothgen(nr_nao, r.data(), orb.data(), rcut_proj, alpha); - - // compare with python script result - std::vector ref = { - 0, - 4.3350439973614511e-05, - 0.00017167638355532129, - 0.00038242839374460959, - 0.0006731078535961104, - 0.0010412661227313883, - 0.0014845037064463464, - 0.0020004694372387291, - 0.0025868596685825239, - 0.0032414174807783905, - 0.0039619318987114734, - 0.0047462371213502306, - 0.0055922117628221264, - 0.006497778104904167, - 0.0074609013607684835, - 0.0084795889498252546, - 0.0095518897835073727, - 0.010675893561843302, - 0.01184973008066669, - 0.013071568549313281 - }; - - for (int i = 0; i < 20; ++i) { - //std::cout< +#include +#include +#include + +namespace +{ + +double radial_norm(const std::vector& r, const std::vector& radial) +{ + std::vector dr(radial.size()); + std::vector integrand(radial.size()); + std::adjacent_difference(r.begin(), r.begin() + radial.size(), dr.begin()); + for (std::size_t i = 0; i < radial.size(); ++i) + { + integrand[i] = r[i] * r[i] * radial[i] * radial[i]; + } + return ModuleBase::Integral::simpson(radial.size(), integrand.data(), &dr[1]); +} + +} // namespace + +TEST(SmoothgenTest, SmoothsAndNormalizesAtCutoff) +{ + const int nr = 121; + const double dr = 0.05; + const double rcut = 4.0; + std::vector r(nr); + std::vector chi(nr); + for (int i = 0; i < nr; ++i) + { + r[i] = i * dr; + chi[i] = std::exp(-0.5 * r[i] * r[i]); + } + + std::vector alpha; + smoothgen(nr, r.data(), chi.data(), rcut, alpha); + + ASSERT_EQ(alpha.size(), 81U); + EXPECT_TRUE(std::all_of(alpha.begin(), alpha.end(), [](const double value) { return std::isfinite(value); })); + EXPECT_NEAR(alpha.back(), 0.0, 1.0e-14); + EXPECT_NEAR(radial_norm(r, alpha), 1.0, 1.0e-10); +} diff --git a/source/source_io/module_hs/unittests/CMakeLists.txt b/source/source_io/module_hs/unittests/CMakeLists.txt index 05bc0102d85..cab6dedb809 100644 --- a/source/source_io/module_hs/unittests/CMakeLists.txt +++ b/source/source_io/module_hs/unittests/CMakeLists.txt @@ -97,6 +97,7 @@ AddTest( ../../../source_basis/module_nao/hydrogen_radials.cpp ../../../source_basis/module_nao/pswfc_radials.cpp ../../../source_basis/module_nao/atomic_radials.cpp + ../../../source_basis/module_nao/smoothgen.cpp ../../../source_basis/module_nao/sphbes_radials.cpp ../../../source_basis/module_nao/two_center_integrator.cpp ../../../source_basis/module_nao/two_center_table.cpp diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index 41ee440a3d7..6cfbecf5089 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -133,7 +133,7 @@ AddTest( if(ENABLE_LCAO) AddTest( TARGET MODULE_IO_to_qo_test - LIBS parameter base device numerical_atomic_orbitals container orb + LIBS parameter base device numerical_atomic_orbitals nao_smoothgen container orb SOURCES to_qo_test.cpp ../module_qo/to_qo_kernel.cpp @@ -173,7 +173,7 @@ add_test(NAME MODULE_IO_read_wf2rho_pw_parallel if(ENABLE_LCAO) AddTest( TARGET MODULE_IO_numerical_basis_test - LIBS parameter base device numerical_atomic_orbitals container orb + LIBS parameter base device numerical_atomic_orbitals nao_smoothgen container orb SOURCES numerical_basis_test.cpp ../module_bessel/numerical_basis_jyjy.cpp ../../source_lcao/center2orb.cpp diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt index b7d0510e0d9..2ee3e1970c2 100644 --- a/source/source_lcao/module_deepks/test/CMakeLists.txt +++ b/source/source_lcao/module_deepks/test/CMakeLists.txt @@ -108,6 +108,7 @@ set(DEEPKS_UNIT_LIBS orb gint numerical_atomic_orbitals + nao_smoothgen symmetry ${TORCH_LIBRARIES} ) diff --git a/source/source_lcao/module_rt/test/CMakeLists.txt b/source/source_lcao/module_rt/test/CMakeLists.txt index 07f0018a68f..02bf3a6a207 100644 --- a/source/source_lcao/module_rt/test/CMakeLists.txt +++ b/source/source_lcao/module_rt/test/CMakeLists.txt @@ -41,7 +41,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_tddft_snap_psibeta_half_test - LIBS parameter base device orb numerical_atomic_orbitals tddft_test_lib + LIBS parameter base device orb numerical_atomic_orbitals nao_smoothgen tddft_test_lib SOURCES snap_psb_half_tddft_test.cpp ../radial_interpolation.cpp ../snap_proj_half_tddft.cpp ../snap_psb_half_tddft.cpp ../../center2orb.cpp ../../center2orb_orb11.cpp diff --git a/source/source_pw/module_proj/CMakeLists.txt b/source/source_pw/module_proj/CMakeLists.txt index bd34f0ddb3f..54282ae771b 100644 --- a/source/source_pw/module_proj/CMakeLists.txt +++ b/source/source_pw/module_proj/CMakeLists.txt @@ -13,6 +13,7 @@ list(APPEND objects radial_proj.cpp + projgen.cpp onsite_proj_print.cpp onsite_proj.cpp onsite_proj_init.cpp diff --git a/source/source_pw/module_proj/onsite_proj_init.cpp b/source/source_pw/module_proj/onsite_proj_init.cpp index 76f0faa7af2..e2c0a4fe411 100644 --- a/source/source_pw/module_proj/onsite_proj_init.cpp +++ b/source/source_pw/module_proj/onsite_proj_init.cpp @@ -9,7 +9,7 @@ #include "source_pw/module_proj/onsite_proj.h" #include "source_pw/module_proj/onsite_proj_print.h" #include "source_pw/module_proj/radial_proj.h" -#include "source_base/projgen.h" +#include "source_basis/module_nao/smoothgen.h" #include "source_base/kernels/math_kernel_op.h" #include "source_base/tool_quit.h" #include "source_base/timer.h" diff --git a/source/source_pw/module_proj/projgen.cpp b/source/source_pw/module_proj/projgen.cpp new file mode 100644 index 00000000000..9fa5de3c7e4 --- /dev/null +++ b/source/source_pw/module_proj/projgen.cpp @@ -0,0 +1,77 @@ +#include "projgen.h" + +#include +#include +#include +#include +#include + +#include "source_base/math_integral.h" +#include "source_base/math_sphbes.h" + +using namespace ModuleBase; + +void projgen(const int l, + const int nr, + const double* r, + const double* chi, + const double rcut, + const int nbes, + std::vector& alpha) +{ + assert(rcut < r[nr - 1]); + assert(std::is_sorted(r, r + nr)); + + std::vector dr(nr); + std::adjacent_difference(r, r + nr, dr.begin()); + + // lower_bound returns the first element that is equal or larger than rcut + int nr_proj = std::distance(r, std::lower_bound(r, r + nr, rcut)) + 1; + + // zeros of spherical Bessel function + std::vector theta(nbes); + Sphbes::sphbes_zeros(l, nbes, theta.data()); + + // z & w vectors (see notes) + std::vector z(nbes); + std::vector w(nbes); + + std::transform(theta.begin(), theta.end(), z.begin(), [rcut, l](double theta_p) { + return 0.5 * std::pow(rcut, 3) * std::pow(Sphbes::sphbesj(l + 1, theta_p), 2); + }); + + // r^2 * chi (independent from p) + std::vector tmp(nr_proj); + std::transform(r, r + nr_proj, chi, tmp.begin(), [](double r_i, double chi_i) { return r_i * r_i * chi_i; }); + + // r^2 * chi * j_l(theta[p] * r / rcut) (dependent on p) + std::vector integrand(nr_proj); + + for (int p = 0; p < nbes; ++p) + { + std::transform(r, r + nr_proj, tmp.begin(), integrand.begin(), [theta, p, rcut, l](double r_i, double tmp_i) { + return tmp_i * Sphbes::sphbesj(l, theta[p] * r_i / rcut); + }); + w[p] = Integral::simpson(nr_proj, integrand.data(), &dr[1]); + } + + // optimal coefficients + std::vector c(nbes, 0.0); + std::transform(w.begin(), w.end(), z.begin(), c.begin(), [](double w_p, double z_p) { return w_p * w_p / z_p; }); + double prefac = 1.0 / std::sqrt(std::accumulate(c.begin(), c.end(), 0.0)); + std::transform(w.begin(), w.end(), z.begin(), c.begin(), [prefac](double w_p, double z_p) { + return prefac * w_p / z_p; + }); + + // new radial function + alpha.resize(nr_proj); + std::fill(alpha.begin(), alpha.end(), 0.0); + for (int i = 0; i < nr_proj; ++i) + { + for (int p = 0; p < nbes; ++p) + { + alpha[i] += c[p] * Sphbes::sphbesj(l, theta[p] * r[i] / rcut); + } + } +} + diff --git a/source/source_base/projgen.h b/source/source_pw/module_proj/projgen.h similarity index 75% rename from source/source_base/projgen.h rename to source/source_pw/module_proj/projgen.h index 5d886ec06c4..99b0f79d917 100644 --- a/source/source_base/projgen.h +++ b/source/source_pw/module_proj/projgen.h @@ -1,5 +1,5 @@ -#ifndef MODULE_NAO_PROJGEN_H_ -#define MODULE_NAO_PROJGEN_H_ +#ifndef MODULE_PROJ_PROJGEN_H +#define MODULE_PROJ_PROJGEN_H #include @@ -25,19 +25,12 @@ * */ void projgen( - const int l, - const int nr, - const double* r, - const double* chi, - const double rcut, - const int nbes, - std::vector& alpha); - -void smoothgen( - const int nr, - const double* r, - const double* chi, - const double rcut, + const int l, + const int nr, + const double* r, + const double* chi, + const double rcut, + const int nbes, std::vector& alpha); #endif \ No newline at end of file diff --git a/source/source_pw/module_proj/unittests/CMakeLists.txt b/source/source_pw/module_proj/unittests/CMakeLists.txt index b9b9154c02b..3c4db87c7ab 100644 --- a/source/source_pw/module_proj/unittests/CMakeLists.txt +++ b/source/source_pw/module_proj/unittests/CMakeLists.txt @@ -3,3 +3,9 @@ AddTest( LIBS parameter base device SOURCES test_radial_proj.cpp ../radial_proj.cpp ) + +AddTest( + TARGET MODULE_PW_projgen_test + LIBS parameter base device + SOURCES test_projgen.cpp ../projgen.cpp +) diff --git a/source/source_base/test/projgen_test.cpp b/source/source_pw/module_proj/unittests/test_projgen.cpp similarity index 64% rename from source/source_base/test/projgen_test.cpp rename to source/source_pw/module_proj/unittests/test_projgen.cpp index 1a70a40cc29..46a096e1861 100644 --- a/source/source_base/test/projgen_test.cpp +++ b/source/source_pw/module_proj/unittests/test_projgen.cpp @@ -1,4 +1,4 @@ -#include "source_base/projgen.h" +#include "source_pw/module_proj/projgen.h" #include "source_base/math_integral.h" @@ -46,24 +46,3 @@ TEST(ProjgenTest, GeneratesNormalizedTruncatedProjector) EXPECT_NEAR(radial_norm(r, alpha), 1.0, 1.0e-10); } -TEST(ProjgenTest, SmoothsAndNormalizesAtCutoff) -{ - const int nr = 121; - const double dr = 0.05; - const double rcut = 4.0; - std::vector r(nr); - std::vector chi(nr); - for (int i = 0; i < nr; ++i) - { - r[i] = i * dr; - chi[i] = std::exp(-0.5 * r[i] * r[i]); - } - - std::vector alpha; - smoothgen(nr, r.data(), chi.data(), rcut, alpha); - - ASSERT_EQ(alpha.size(), 81U); - EXPECT_TRUE(std::all_of(alpha.begin(), alpha.end(), [](const double value) { return std::isfinite(value); })); - EXPECT_NEAR(alpha.back(), 0.0, 1.0e-14); - EXPECT_NEAR(radial_norm(r, alpha), 1.0, 1.0e-10); -}