From 42be6aff3f702d05882fd12681d62d2e822d8a87 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 16:56:38 +0800 Subject: [PATCH 01/14] cherry pick from TT & LRT's commit for socket mode --- doc/.LaTeX/Socket.tex | 215 ++++++ doc/.LaTeX/mlff/Manual_mlff.tex | 73 ++ src/electronicGroundState.c | 35 +- src/main.c | 5 +- src/main_socket.c | 80 +++ src/md.c | 298 +++++++- src/relax.c | 34 +- src/socket/conf.h | 9 + src/socket/driver.c | 1164 +++++++++++++++++++++++++++++++ src/socket/driver.h | 74 ++ src/socket/libinetsocket.c | 1106 +++++++++++++++++++++++++++++ src/socket/libinetsocket.h | 97 +++ src/socket/libunixsocket.c | 515 ++++++++++++++ src/socket/libunixsocket.h | 73 ++ 14 files changed, 3729 insertions(+), 49 deletions(-) create mode 100644 doc/.LaTeX/Socket.tex create mode 100644 doc/.LaTeX/mlff/Manual_mlff.tex create mode 100644 src/main_socket.c create mode 100644 src/socket/conf.h create mode 100644 src/socket/driver.c create mode 100644 src/socket/driver.h create mode 100644 src/socket/libinetsocket.c create mode 100644 src/socket/libinetsocket.h create mode 100644 src/socket/libunixsocket.c create mode 100644 src/socket/libunixsocket.h diff --git a/doc/.LaTeX/Socket.tex b/doc/.LaTeX/Socket.tex new file mode 100644 index 00000000..a8635b53 --- /dev/null +++ b/doc/.LaTeX/Socket.tex @@ -0,0 +1,215 @@ + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks,c]{} \label{Socket} + +\begin{center} +\Huge \textbf{Socket communication in SPARC} +\end{center} + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\textbf{Brief introduction}} \label{SOCKET_INTRO} + + The socket communication layer allows SPARC to be controlled by a socket server that is compatible with + \href{https://ipi-code.org/i-pi/introduction.html}{i-PI} protocol. The SPARC source code should be compiled with the + \texttt{USE\_SOCKET=1} option. + + To start a SPARC program with socket interface, use either option: + \begin{itemize} + \item Option 1: Specify socket address in command line +\begin{itemize} + \item To start an INET socket, use\\ + \texttt{\$ mpirun -n 8 ./lib/sparc -socket localhost:12345 -name filename} + + \item To start a UNIX socket, use\\ + \texttt{\$ mpirun -n 8 ./lib/sparc -socket /tmp/sparc.socket:unix -name filename} + \end{itemize} + + \item Option 2: Provide \texttt{.inpt} parameters + + Visit the following pages for parameter specifications. + Any parameter \texttt{SPARC\_*} in the \texttt{.inpt} file + will overwrite the settings from command line. + \end{itemize} + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\textbf{Output files}} \label{SOCKET_OUTPUT} + + The following files will be generated in the socket mode and may slightly + differ from standard SPARC outputs: + \begin{itemize} + \item \texttt{.out}: All the SCF steps will be written to \texttt{.out} file, similar to a relaxation / MD calculation. + \item \texttt{.static}: It is concatenated from all single point + steps. In addition to normal \texttt{.static} output, the lattice + information are also recorded. + \end{itemize} + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\texttt{SOCKET\_FLAG}} \label{SOCKET_FLAG} +\vspace*{-12pt} +\begin{columns} +\column{0.4\linewidth} +\begin{block}{Type} +Integer +\end{block} + +\begin{block}{Default} +0 +\end{block} + +\column{0.4\linewidth} +\begin{block}{Unit} +No unit +\end{block} + +\begin{block}{Example} +\texttt{SOCKET\_FLAG}: 1 +\end{block} +\end{columns} + +\begin{block}{Description} + Flag for starting the socket communication layer. It is equivalent to the \texttt{-socket} switch of command line. + + Setting \texttt{SOCKET\_FLAG: 1} will disable \texttt{MD\_FLAG} and \texttt{RELAX\_FLAG}. + +\end{block} + + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\texttt{SOCKET\_HOST}} \label{SOCKET_HOST} +\vspace*{-12pt} +\begin{columns} +\column{0.4\linewidth} +\begin{block}{Type} +String +\end{block} + +\begin{block}{Default} +localhost +\end{block} + +\column{0.4\linewidth} +\begin{block}{Unit} +No unit +\end{block} + +\begin{block}{Example} +\texttt{SOCKET\_HOST}: 127.0.0.1 +\end{block} +\end{columns} + +\begin{block}{Description} + Host name of the socket address that SPARC listens to. If it's an + INET socket, it is the address of the interface. For a UNIX socket, + it is the filename of the socket file + (e.g. \texttt{/tmp/sparc.socket}). +\end{block} + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\texttt{SOCKET\_PORT}} \label{SOCKET_PORT} +\vspace*{-12pt} +\begin{columns} +\column{0.4\linewidth} +\begin{block}{Type} +Integer +\end{block} + +\begin{block}{Default} +None +\end{block} + +\column{0.4\linewidth} +\begin{block}{Unit} +No unit +\end{block} + +\begin{block}{Example} +\texttt{SOCKET\_PORT}: 12345 +\end{block} +\end{columns} + +\begin{block}{Description} + When SPARC connects to an INET socket server, it is the port number. The \texttt{SOCKET\_PORT} has no effect for a UNIX socket. +\end{block} + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\texttt{SOCKET\_INET}} \label{SOCKET_INET} +\vspace*{-12pt} +\begin{columns} +\column{0.4\linewidth} +\begin{block}{Type} +Integer +\end{block} + +\begin{block}{Default} +None +\end{block} + +\column{0.4\linewidth} +\begin{block}{Unit} +No unit +\end{block} + +\begin{block}{Example} +\texttt{SOCKET\_INET}: 0 +\end{block} +\end{columns} + +\begin{block}{Description} + 1 for INET socket, 0 for UNIX socket. If no set in \texttt{.inpt} + file, its value is determined by the \texttt{-socket address:port} + command line switch. +\end{block} + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\begin{frame}[allowframebreaks]{\texttt{SOCKET\_MAX\_NITER}} \label{SOCKET_MAX_NITER} +\vspace*{-12pt} +\begin{columns} +\column{0.4\linewidth} +\begin{block}{Type} +Integer +\end{block} + +\begin{block}{Default} +10000 +\end{block} + +\column{0.4\linewidth} +\begin{block}{Unit} +No unit +\end{block} + +\begin{block}{Example} +\texttt{SOCKET\_MAX\_NITER}: 10000 +\end{block} +\end{columns} + +\begin{block}{Description} + Maximum number of ionic SCF steps in the socket mode. As a socket + client, SPARC will terminate after \texttt{SOCKET\_MAX\_NITER} steps + are called. +\end{block} + + +\end{frame} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \ No newline at end of file diff --git a/doc/.LaTeX/mlff/Manual_mlff.tex b/doc/.LaTeX/mlff/Manual_mlff.tex new file mode 100644 index 00000000..ed0a40e7 --- /dev/null +++ b/doc/.LaTeX/mlff/Manual_mlff.tex @@ -0,0 +1,73 @@ +\documentclass[xcolor=dvipsnames,t]{beamer} +%\usepackage[utf8]{inputenc} +\usepackage{verbatim} % for comment +\usetheme{Madrid} +\usecolortheme{seahorse} +\usepackage{beamerouterthemesplit} +\usepackage[T1]{fontenc} + +\setbeamerfont{institute}{size=\normalsize} + +\title[] {SPARC - MLFF \\ {\normalsize Machine-learned force fields} \\ {\small User guide}} +%\author{} +%\author{Qimen Xu, Abhiraj Sharma, Phanish Suryanarayana, } + +\institute[] +{ +Material Physics \& Mechanics Group \\ +PI: Phanish Suryanarayana \\ +%Main Developers: Shashikant Kumar, Abhiraj Sharma\\ +%Collaborators: J.E. Pask (LLNL)\\ +Georgia Institute of Technology \\ +\hyperlink{Contributors}{Contributors}\\ +\hyperlink{Citation}{Citation} \\ +\hyperlink{Acknowledgements}{Acknowledgements} +} + +\date{} +\setbeamertemplate{frametitle continuation}{} +\setbeamertemplate{navigation symbols}{} +\setbeamercolor{block title}{bg=Apricot!50,fg=black} +%\addtobeamertemplate{block begin}{\vskip -\smallskipamount}{} +%\addtobeamertemplate{block end}{}{\vskip -\smallskipamount} +\addtobeamertemplate{block begin}{\vspace*{-0.6pt}}{} +\addtobeamertemplate{block end}{}{\vspace*{-0.6pt}} +\hypersetup{colorlinks, +citecolor=violet, +linkcolor=blue, +menucolor=white, +anchorcolor=yellow +filecolor=pink, +} + + +\setbeamertemplate{footline}{ + \hbox{% + \begin{beamercolorbox}[wd=\paperwidth,ht=1ex,dp=1.5ex,leftskip=2ex,rightskip=2ex]{page footer}% + \usebeamerfont{title in head/foot}% + \insertshorttitle \hfill + \insertsection \hfill + \insertframenumber{} / \inserttotalframenumber + \end{beamercolorbox}}% +} + +\let\otp\titlepage +\renewcommand{\titlepage}{\otp\addtocounter{framenumber}{-1}} + + +%\includeonly{Introduction,Stress,MD,Optimization,Print} +%\includeonly{Introduction} + +\begin{document} + +%\frame{\titlepage} + +\begin{frame}[plain] + \titlepage +\end{frame} + +\include{Introduction} +\include{mlff} + + +\end{document} diff --git a/src/electronicGroundState.c b/src/electronicGroundState.c index bf013af6..57c289ac 100644 --- a/src/electronicGroundState.c +++ b/src/electronicGroundState.c @@ -53,6 +53,7 @@ #include "sqNlocVecRoutines.h" #include "ofdft.h" #include "printing.h" +#include "sparc_mlff_interface.h" #ifdef USE_EVA_MODULE #include "ExtVecAccel/ExtVecAccel.h" @@ -62,6 +63,16 @@ #define min(x,y) ((x)<(y)?(x):(y)) #define TEMP_TOL 1e-12 +/** + * @brief Calculate properties for fixed atom positions using DFT or MLFF + */ +void Calculate_Properties(SPARC_OBJ *pSPARC) { + if (pSPARC->mlff_flag < 1) + Calculate_electronicGroundState(pSPARC); + else + MLFF_call(pSPARC); +} + /** * @brief Calculate the ground state energy and forces for fixed atom positions. */ @@ -86,6 +97,10 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) { SCF_ind = pSPARC->RelaxCount + pSPARC->restartCount + (pSPARC->RestartFlag == 0); else SCF_ind = 1; +#ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) + SCF_ind = pSPARC->SocketSCFCount; +#endif if (pSPARC->REFERENCE_CUTOFF > 0.5*nn) { printf("\nWARNING: REFERENCE_CUFOFF (%.6f Bohr) > 1/2 nn (nearest neighbor) distance (%.6f Bohr) in SCF#%d\n", pSPARC->REFERENCE_CUTOFF, 0.5*nn, SCF_ind); @@ -144,6 +159,8 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) { } fclose(output_fp); // for static calculation, print energy to .static file + // if socket mode is activated, each static scf step is + // separated by comment lines if (pSPARC->MDFlag == 0 && pSPARC->RelaxFlag == 0) { if (pSPARC->PrintForceFlag == 1 || pSPARC->PrintAtomPosFlag == 1) { static_fp = fopen(pSPARC->StaticFilename,"a"); @@ -164,7 +181,7 @@ void Calculate_electronicGroundState(SPARC_OBJ *pSPARC) { // calculate atom magnetization if (pSPARC->spin_typ) CalculateAtomMag(pSPARC); - + // write forces into .static file if required if (rank == 0 && pSPARC->Verbosity && pSPARC->PrintForceFlag == 1 && pSPARC->MDFlag == 0 && pSPARC->RelaxFlag == 0) { static_fp = fopen(pSPARC->StaticFilename,"a"); @@ -553,8 +570,13 @@ void scf(SPARC_OBJ *pSPARC) else if(pSPARC->RelaxFlag >= 1) fprintf(output_fp," Self Consistent Field (SCF#%d) \n", pSPARC->RelaxCount + pSPARC->restartCount + (pSPARC->RestartFlag == 0)); +#ifdef USE_SOCKET + else if(pSPARC->SocketFlag == 1) + fprintf(output_fp," Self Consistent Field (SCF#%d) \n", pSPARC->SocketSCFCount); +#endif else - fprintf(output_fp," Self Consistent Field (SCF#%d) \n",1); + fprintf(output_fp," Self Consistent Field (SCF#%d) \n",1); + if(pSPARC->spin_typ == 0) fprintf(output_fp,"===================================================================\n"); @@ -614,8 +636,8 @@ void scf(SPARC_OBJ *pSPARC) if (rank == 0) { printf("rank = %d, Veff calculation and Bcast (non-blocking) took %.3f ms\n",rank,(t2-t1)*1e3); } - #endif - + #endif + scf_loop(pSPARC); if (pSPARC->usefock > 0) { pSPARC->usefock ++; @@ -738,7 +760,8 @@ void scf_loop(SPARC_OBJ *pSPARC) { // calculate xc potential (LDA, PW92), "Vxc" Calculate_Vxc(pSPARC); - pSPARC->countPotentialCalculate++; + pSPARC->countPotentialCalculate++; + #ifdef DEBUG t2 = MPI_Wtime(); if(!rank) printf("rank = %d, Calculating Vxc took %.3f ms\n", rank, (t2 - t1) * 1e3); @@ -1398,4 +1421,4 @@ void CalculateAtomMag(SPARC_OBJ *pSPARC) } MPI_Allreduce(MPI_IN_PLACE, pSPARC->AtomMag, pSPARC->n_atom*ncol, MPI_DOUBLE, MPI_SUM, pSPARC->dmcomm_phi); -} \ No newline at end of file +} diff --git a/src/main.c b/src/main.c index c847d6a6..0eb9a17b 100644 --- a/src/main.c +++ b/src/main.c @@ -47,9 +47,8 @@ int main(int argc, char *argv[]) { else if (SPARC.RelaxFlag != 0) main_Relax(&SPARC); else - Calculate_electronicGroundState(&SPARC); - - + Calculate_Properties(&SPARC); + // Calculate_electronicGroundState(&SPARC); Finalize(&SPARC); diff --git a/src/main_socket.c b/src/main_socket.c new file mode 100644 index 00000000..59e733aa --- /dev/null +++ b/src/main_socket.c @@ -0,0 +1,80 @@ +/** + * @file main_socket.c + * @brief This file contains the main function for SPARC + * + * @authors Qimen Xu + * Tian Tian + * Abhiraj Sharma + * Phanish Suryanarayana + * + * Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech. + */ + +#include +#include +#include +#include + +#include "initialization.h" +#include "md.h" +#include "relax.h" +//#include "tools.h" +#include "finalization.h" +#include "isddft.h" +#include "electronicGroundState.h" + +#ifdef USE_SOCKET +#include "driver.h" +#endif + + +int main(int argc, char *argv[]) { + // set up MPI + MPI_Init(&argc, &argv); + // get communicator size and my rank + MPI_Comm comm = MPI_COMM_WORLD; + int nproc, rank; + MPI_Comm_size(comm, &nproc); + MPI_Comm_rank(comm, &rank); + + SPARC_OBJ SPARC; + + double t1, t2; + + MPI_Barrier(MPI_COMM_WORLD); + // start timer + t1 = MPI_Wtime(); SPARC.time_start = t1; + + // Read files and initialize + Initialize(&SPARC, argc, argv); + + if (SPARC.MDFlag == 1) + main_MD(&SPARC); + else if (SPARC.RelaxFlag != 0) + main_Relax(&SPARC); + else +#ifdef USE_SOCKET + if (SPARC.SocketFlag == 1) + main_Socket(&SPARC); + else + Calculate_Properties(&SPARC); + //Calculate_electronicGroundState(&SPARC); +#else + Calculate_Properties(&SPARC); + //Calculate_electronicGroundState(&SPARC); +#endif + Finalize(&SPARC); + + MPI_Barrier(MPI_COMM_WORLD); + // end timer + t2 = MPI_Wtime(); + if (rank == 0) { + printf("The program took %.3f s.\n", t2 - t1); + } + // ensure stdout flushed to prevent Finalize hang + fflush(stdout); + + // finalize MPI + MPI_Finalize(); + return 0; +} diff --git a/src/md.c b/src/md.c index 42eda887..7f9bd2ef 100644 --- a/src/md.c +++ b/src/md.c @@ -72,15 +72,9 @@ void main_MD(SPARC_OBJ *pSPARC) { } } - - - // Initialize the MD and MLFF for the very first step only - if (pSPARC->mlff_flag==0){ - Calculate_electronicGroundState(pSPARC); - Initialize_MD(pSPARC); - } else { - MLFF_main(pSPARC); - } + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); + Initialize_MD(pSPARC); pSPARC->MD_maxStep = pSPARC->restartCount + pSPARC->MD_Nstep; @@ -503,12 +497,9 @@ void NVT_NH(SPARC_OBJ *pSPARC) { elecDensExtrapolation(pSPARC); // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); - // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF - if (pSPARC->mlff_flag == 0){ - Calculate_electronicGroundState(pSPARC); - } else { - MLFF_call(pSPARC); - } + // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; // Second step velocity Verlet VVerlet2(pSPARC); @@ -648,12 +639,9 @@ void NVE(SPARC_OBJ *pSPARC) { elecDensExtrapolation(pSPARC); // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); - // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF - if (pSPARC->mlff_flag == 0){ - Calculate_electronicGroundState(pSPARC); - } else { - MLFF_call(pSPARC); - } + // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; // Leapfrog step (part-2) Leapfrog_part2(pSPARC); @@ -727,13 +715,9 @@ void NVK_G(SPARC_OBJ *pSPARC) { // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); - - // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF - if (pSPARC->mlff_flag == 0){ - Calculate_electronicGroundState(pSPARC); - } else { - MLFF_call(pSPARC); - } + // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPSPARC); pSPARC->elecgs_Count++; // Calculate velocity at next full time step @@ -870,14 +854,20 @@ void NPT_NH (SPARC_OBJ *pSPARC) { MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&pSPARC->pres_i, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); +<<<<<<< variant A +>>>>>>> variant B +======= end if ((pSPARC->MDCount != 1) || (pSPARC->RestartFlag == 1)) { // Update velocity of particles in the second half timestep AccelVelocityParticle(pSPARC); // Update velocity of virtual thermo and baro variables in the second half timestep IsoPress(pSPARC); } +<<<<<<< variant A +>>>>>>> variant B +======= end // Update velocity of virtual thermo and baro variables in the first half timestep IsoPress(pSPARC); // Update velocity of particles in the first half timestep @@ -886,17 +876,26 @@ void NPT_NH (SPARC_OBJ *pSPARC) { PositionParticleCell(pSPARC); // Reinitialize mesh size and related variables after changing size of cell reinitialize_mesh_NPT(pSPARC); +<<<<<<< variant A +>>>>>>> variant B +======= end // Charge extrapolation (for better rho_guess) elecDensExtrapolation(pSPARC); // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); +<<<<<<< variant A // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF if (pSPARC->mlff_flag == 0){ Calculate_electronicGroundState(pSPARC); } else { MLFF_call(pSPARC); } +>>>>>>> variant B + // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); +======= end pSPARC->elecgs_Count++; #ifdef DEBUG // Calculate Hamiltonian of the system. @@ -905,7 +904,10 @@ void NPT_NH (SPARC_OBJ *pSPARC) { printf("\nend NPT timestep %d\n", pSPARC->MDCount + 1); } #endif +<<<<<<< variant A +>>>>>>> variant B +======= end } /* @@ -1137,7 +1139,11 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) double t1, t2; #endif +<<<<<<< variant A int p, i; +>>>>>>> variant B + int p; +======= end // scaling factor double scal = pSPARC->scale; // the ratio between length #ifdef DEBUG @@ -1150,9 +1156,13 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->delta_y = pSPARC->range_y/(pSPARC->numIntervals_y); pSPARC->delta_z = pSPARC->range_z/(pSPARC->numIntervals_z); +<<<<<<< variant A pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z * pSPARC->Jacbdet; +>>>>>>> variant B +======= end +<<<<<<< variant A #ifdef DEBUG if (!rank) { // printf("Volume: %12.6f\n", vol); @@ -1163,8 +1173,18 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) if (i%3==2 && i>0) printf("\n"); } printf("\n"); +>>>>>>> variant B + if (pSPARC->cell_typ <= 20) { + pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z * pSPARC->Jacbdet; + } else if (pSPARC->cell_typ > 20 && pSPARC->cell_typ <= 30) { + pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z; + // pSPARC->twist = pSPARC->twistpercell/pSPARC->range_z; +======= end } +<<<<<<< variant A #endif +>>>>>>> variant B +======= end int FDn = pSPARC->order / 2; @@ -1179,6 +1199,7 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->D1_stencil_coeffs_z[p] = pSPARC->FDweights_D1[p] * dz_inv; } +<<<<<<< variant A // 2nd derivative weights including mesh double dx2_inv, dy2_inv, dz2_inv; dx2_inv = 1.0 / (pSPARC->delta_x * pSPARC->delta_x); @@ -1208,7 +1229,21 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->D1_stencil_coeffs_zy[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // d/dz(2*T_23 df/dy) used in d/dz(2*T_12 df/dx + 2*T_23 df/dy) } } - +>>>>>>> variant B + if (pSPARC->CyclixFlag) { + LapStencil_cyclix(pSPARC); + Integration_weights_cyclix(pSPARC, pSPARC->Intgwt_kpttopo, pSPARC->DMVertices_kptcomm[0], pSPARC->Nx_d_kptcomm, pSPARC->Ny_d_kptcomm, pSPARC->Nz_d_kptcomm); + Integration_weights_cyclix(pSPARC, pSPARC->Intgwt_psi, pSPARC->DMVertices_dmcomm[0], pSPARC->Nx_d_dmcomm, pSPARC->Ny_d_dmcomm, pSPARC->Nz_d_dmcomm); + Integration_weights_cyclix(pSPARC, pSPARC->Intgwt_phi, pSPARC->DMVertices[0], pSPARC->Nx_d, pSPARC->Ny_d, pSPARC->Nz_d); + } else { + // 2nd derivative weights including mesh + double dx2_inv, dy2_inv, dz2_inv; + dx2_inv = 1.0 / (pSPARC->delta_x * pSPARC->delta_x); + dy2_inv = 1.0 / (pSPARC->delta_y * pSPARC->delta_y); + dz2_inv = 1.0 / (pSPARC->delta_z * pSPARC->delta_z); +======= end + +<<<<<<< variant A // maximum eigenvalue of -0.5 * Lap (only with periodic boundary conditions) if(pSPARC->cell_typ == 0) { #ifdef DEBUG @@ -1225,7 +1260,32 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->MaxEigVal_mhalfLap += 2.0 * (pSPARC->D2_stencil_coeffs_x[p] * cos(M_PI*p*scal_x) + pSPARC->D2_stencil_coeffs_y[p] * cos(M_PI*p*scal_y) + pSPARC->D2_stencil_coeffs_z[p] * cos(M_PI*p*scal_z)); +>>>>>>> variant B + // Stencil coefficients for mixed derivatives + if (pSPARC->cell_typ == 0) { + for (p = 0; p < FDn + 1; p++) { + pSPARC->D2_stencil_coeffs_x[p] = pSPARC->FDweights_D2[p] * dx2_inv; + pSPARC->D2_stencil_coeffs_y[p] = pSPARC->FDweights_D2[p] * dy2_inv; + pSPARC->D2_stencil_coeffs_z[p] = pSPARC->FDweights_D2[p] * dz2_inv; + } + } else if (pSPARC->cell_typ > 10 && pSPARC->cell_typ < 20) { + for (p = 0; p < FDn + 1; p++) { + pSPARC->D2_stencil_coeffs_x[p] = pSPARC->lapcT[0] * pSPARC->FDweights_D2[p] * dx2_inv; + pSPARC->D2_stencil_coeffs_y[p] = pSPARC->lapcT[4] * pSPARC->FDweights_D2[p] * dy2_inv; + pSPARC->D2_stencil_coeffs_z[p] = pSPARC->lapcT[8] * pSPARC->FDweights_D2[p] * dz2_inv; + pSPARC->D2_stencil_coeffs_xy[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dx_inv; // 2*T_12 d/dx(df/dy) + pSPARC->D2_stencil_coeffs_xz[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dx_inv; // 2*T_13 d/dx(df/dz) + pSPARC->D2_stencil_coeffs_yz[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // 2*T_23 d/dy(df/dz) + pSPARC->D1_stencil_coeffs_xy[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dy_inv; // d/dx(2*T_12 df/dy) used in d/dx(2*T_12 df/dy + 2*T_13 df/dz) + pSPARC->D1_stencil_coeffs_yx[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dx_inv; // d/dy(2*T_12 df/dx) used in d/dy(2*T_12 df/dx + 2*T_23 df/dz) + pSPARC->D1_stencil_coeffs_xz[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dz_inv; // d/dx(2*T_13 df/dz) used in d/dx(2*T_12 df/dy + 2*T_13 df/dz) + pSPARC->D1_stencil_coeffs_zx[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dx_inv; // d/dz(2*T_13 df/dx) used in d/dz(2*T_13 df/dz + 2*T_23 df/dy) + pSPARC->D1_stencil_coeffs_yz[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dz_inv; // d/dy(2*T_23 df/dz) used in d/dy(2*T_12 df/dx + 2*T_23 df/dz) + pSPARC->D1_stencil_coeffs_zy[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // d/dz(2*T_23 df/dy) used in d/dz(2*T_12 df/dx + 2*T_23 df/dy) + } +======= end } +<<<<<<< variant A pSPARC->MaxEigVal_mhalfLap *= -0.5; #ifdef DEBUG t2 = MPI_Wtime(); @@ -1233,7 +1293,10 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->MaxEigVal_mhalfLap, (t2-t1)*1e3); #endif } +>>>>>>> variant B +======= end +<<<<<<< variant A double h_eff = 0.0; if (fabs(pSPARC->delta_x - pSPARC->delta_y) < 1E-12 && fabs(pSPARC->delta_y - pSPARC->delta_z) < 1E-12) { @@ -1242,14 +1305,51 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) // find effective mesh s.t. it has same spectral width h_eff = sqrt(3.0 / (dx2_inv + dy2_inv + dz2_inv)); } +>>>>>>> variant B + // maximum eigenvalue of -0.5 * Lap (only with periodic boundary conditions) + if(pSPARC->cell_typ == 0) { + #ifdef DEBUG + t1 = MPI_Wtime(); + #endif + pSPARC->MaxEigVal_mhalfLap = pSPARC->D2_stencil_coeffs_x[0] + + pSPARC->D2_stencil_coeffs_y[0] + + pSPARC->D2_stencil_coeffs_z[0]; + double scal_x, scal_y, scal_z; + scal_x = (pSPARC->Nx - pSPARC->Nx % 2) / (double) pSPARC->Nx; + scal_y = (pSPARC->Ny - pSPARC->Ny % 2) / (double) pSPARC->Ny; + scal_z = (pSPARC->Nz - pSPARC->Nz % 2) / (double) pSPARC->Nz; + for (int p = 1; p < FDn + 1; p++) { + pSPARC->MaxEigVal_mhalfLap += 2.0 * (pSPARC->D2_stencil_coeffs_x[p] * cos(M_PI*p*scal_x) + + pSPARC->D2_stencil_coeffs_y[p] * cos(M_PI*p*scal_y) + + pSPARC->D2_stencil_coeffs_z[p] * cos(M_PI*p*scal_z)); + } + pSPARC->MaxEigVal_mhalfLap *= -0.5; + #ifdef DEBUG + t2 = MPI_Wtime(); + if (!rank) printf("Max eigenvalue of -0.5*Lap is %.13f, time taken: %.3f ms\n", + pSPARC->MaxEigVal_mhalfLap, (t2-t1)*1e3); + #endif + } +======= end +<<<<<<< variant A // find Chebyshev polynomial degree based on max eigenvalue (spectral width) if (pSPARC->ChebDegree < 0) { pSPARC->ChebDegree = Mesh2ChebDegree(h_eff); #ifdef DEBUG if (!rank && h_eff < 0.1) { printf("WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n"); +>>>>>>> variant B + double h_eff = 0.0; + if (fabs(pSPARC->delta_x - pSPARC->delta_y) < 1E-12 && + fabs(pSPARC->delta_y - pSPARC->delta_z) < 1E-12) { + h_eff = pSPARC->delta_x; + } else { + // find effective mesh s.t. it has same spectral width + h_eff = sqrt(3.0 / (dx2_inv + dy2_inv + dz2_inv)); +======= end } +<<<<<<< variant A if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree); #endif } else { @@ -1257,18 +1357,49 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) if (!rank) printf("Chebyshev polynomial degree (provided by user): npl = %d\n",pSPARC->ChebDegree); #endif } +>>>>>>> variant B +======= end +<<<<<<< variant A // default Kerker tolerance if (pSPARC->TOL_PRECOND < 0.0) { // kerker tol not provided by user pSPARC->TOL_PRECOND = (h_eff * h_eff) * 1e-3; } +>>>>>>> variant B + // find Chebyshev polynomial degree based on max eigenvalue (spectral width) + if (pSPARC->ChebDegree < 0) { + pSPARC->ChebDegree = Mesh2ChebDegree(h_eff); + #ifdef DEBUG + if (!rank && h_eff < 0.1) { + printf("#WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enough!\n"); + } + if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree); + #endif + } else { + #ifdef DEBUG + if (!rank) printf("Chebyshev polynomial degree (provided by user): npl = %d\n",pSPARC->ChebDegree); + #endif + } +======= end +<<<<<<< variant A +>>>>>>> variant B + // default Kerker tolerance + if (pSPARC->TOL_PRECOND < 0.0) { // kerker tol not provided by user + pSPARC->TOL_PRECOND = (h_eff * h_eff) * 1e-3; + } + } +======= end // re-calculate k-point grid Calculate_kpoints(pSPARC); // re-calculate local k-points array +<<<<<<< variant A if (pSPARC->Nkpts >= 1 && pSPARC->kptcomm_index != -1) { +>>>>>>> variant B + if (pSPARC->Nkpts >= 1 && pSPARC->kptcomm_index != -1 && (!pSPARC->sqAmbientFlag)) { +======= end Calculate_local_kpoints(pSPARC); } @@ -1282,12 +1413,26 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) if (rank == 0) printf("Calculating rb for all atom types took %.3f ms\n",(t2-t1)*1000); #endif +<<<<<<< variant A +>>>>>>> variant B + // Rescale planewaves for spDFT + if (pSPARC->spDFT_Flag == 1) { + spDFT_planewaves(pSPARC); + } +======= end + +<<<<<<< variant A +>>>>>>> variant B +======= end // write reinitialized parameters into output file if (rank == 0) { write_output_reinit_NPT(pSPARC); } +<<<<<<< variant A +>>>>>>> variant B +======= end } @@ -1355,17 +1500,26 @@ void NPT_NP (SPARC_OBJ *pSPARC) { updatePosition(pSPARC); // Reinitialize mesh size and related variables after changing size of cell reinitialize_mesh_NPT(pSPARC); +<<<<<<< variant A +>>>>>>> variant B +======= end // Charge extrapolation (for better rho_guess) elecDensExtrapolation(pSPARC); // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); +<<<<<<< variant A // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF if (pSPARC->mlff_flag == 0){ Calculate_electronicGroundState(pSPARC); } else { MLFF_call(pSPARC); } +>>>>>>> variant B + // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); +======= end pSPARC->elecgs_Count++; #ifdef DEBUG if (!rank) printf("\nend NPT_NP timestep %d\n", pSPARC->MDCount + 1); @@ -1483,6 +1637,12 @@ void updateMomentum_FirstHalf(SPARC_OBJ *pSPARC) { innerControlStress[1] = innerControlStress[0]; innerControlStress[2] = innerControlStress[0]; } +<<<<<<< variant A +>>>>>>> variant B + if (pSPARC->CyclixFlag) { // cyclix on z, unit Ha/Bohr. It is necessary to convert it to Ha/Bohr^3. + innerControlStress[2] /= M_PI * (pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0)) ; + } +======= end // update momentum of barostat variables in a half step for (int i = 0; i < 3; i++){ if (pSPARC->NPTscaleVecs[i] == 1) { @@ -1493,7 +1653,11 @@ void updateMomentum_FirstHalf(SPARC_OBJ *pSPARC) { pSPARC->Pm_NPT_NP[i] -= pSPARC->MD_dt * pSPARC->S_NPT_NP / 2.0 * PmA[i]; #ifdef DEBUG if (rank == 0){ +<<<<<<< variant A // printf("pSPARC->pres is %12.9f, pSPARC->pres_i is %12.9f, pSPARC->volumeCell is %12.9f, pSPARC->G_NPT_NP[%d] is %12.9f\n", pSPARC->pres, pSPARC->pres_i, pSPARC->volumeCell, i, pSPARC->G_NPT_NP[i]); +>>>>>>> variant B + printf("innerControlStress[%d] is %12.9f, pSPARC->volumeCell is %12.9f, pSPARC->G_NPT_NP[%d] is %12.9f\n", i, innerControlStress[i], pSPARC->volumeCell, i, pSPARC->G_NPT_NP[i]); +======= end printf("PmA[%d] is %12.9f\n", i, PmA[i]); printf("pSPARC->Pm_NPT_NP[%d] in 1st half step is %12.9f\n", i, pSPARC->Pm_NPT_NP[i]); } @@ -1557,7 +1721,13 @@ void updateMomentum_SecondHalf(SPARC_OBJ *pSPARC) { innerControlStress[1] = innerControlStress[0]; innerControlStress[2] = innerControlStress[0]; } +<<<<<<< variant A +>>>>>>> variant B + if (pSPARC->CyclixFlag) { // cyclix on z, unit Ha/Bohr. It is necessary to convert it to Ha/Bohr^3. + innerControlStress[2] /= M_PI * (pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0)) ; + } +======= end while (judge == 0) { timeIter++; KbaroTmp = 0; @@ -1721,7 +1891,17 @@ void updatePosition(SPARC_OBJ *pSPARC) { pSPARC->range_y = sqrt(pSPARC->G_NPT_NP[1]); double scalez = sqrt(pSPARC->G_NPT_NP[2]) / pSPARC->range_z; pSPARC->range_z = sqrt(pSPARC->G_NPT_NP[2]); +<<<<<<< variant A pSPARC->volumeCell = pSPARC->Jacbdet*pSPARC->range_x*pSPARC->range_y*pSPARC->range_z; +>>>>>>> variant B + if (pSPARC->CyclixFlag) { + pSPARC->twist = pSPARC->twistpercell / pSPARC->range_z; + pSPARC->volumeCell = pSPARC->range_x * ((pSPARC->xin + pSPARC->xout)/2.0) * pSPARC->range_y * pSPARC->range_z; + } + else { + pSPARC->volumeCell = pSPARC->Jacbdet*pSPARC->range_x*pSPARC->range_y*pSPARC->range_z; + } +======= end if (pSPARC->NPTscaleVecs[0] == 1) pSPARC->range_x_velo = G3[0] * Stemp / (2.0*pSPARC->NPT_NP_bmass*pow(pSPARC->volumeCell, 2.0)) / pSPARC->range_x; else @@ -1768,7 +1948,11 @@ void Cart2nonCart(double *gradT, double *carCoord, double *nonCarCoord) { } /* +<<<<<<< variant A * @ brief: function to check if the atoms are too close to the boundary in case of bounded domain or to each other in general +>>>>>>> variant B + @ brief: function to wrap around atom positions that lie outside main domain in case of PBC and check if the atoms are too close to the boundary in case of bounded domain +======= end */ void Check_atomlocation(SPARC_OBJ *pSPARC) { int rank, ityp, i, atm, atm2, count, dir = 0, maxdir = 3, BC; @@ -1816,9 +2000,17 @@ void Check_atomlocation(SPARC_OBJ *pSPARC) { wraparound_dynamics(pSPARC, pSPARC->atom_pos, 1); } +<<<<<<< variant A /* * @ brief: function to wraparound atom positions for PBC */ +>>>>>>> variant B + + + + + +======= end void wraparound_dynamics(SPARC_OBJ *pSPARC, double *coord, int opt) { int rank, atm, dir = 0, maxdir = 3, BC; @@ -1873,9 +2065,13 @@ void wraparound_dynamics(SPARC_OBJ *pSPARC, double *coord, int opt) { } } +<<<<<<< variant A /* * @ brief: function to wraparound velocities in MD and displacement vectors in relaxation for PBC */ +>>>>>>> variant B + +======= end void wraparound_velocity(SPARC_OBJ *pSPARC, double shift, int dir, int loc) { if (dir == 1){ @@ -1904,6 +2100,12 @@ void wraparound_velocity(SPARC_OBJ *pSPARC, double shift, int dir, int loc) { } +<<<<<<< variant A +>>>>>>> variant B + + + +======= end /* @ brief: function to write all relevant DFT quantities generated during MD simulation */ @@ -2153,7 +2355,11 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) { } // Average and minimum distance +<<<<<<< variant A //*avgdis = pow((pSPARC->range_x * pSPARC->range_y * pSPARC->range_z)/pSPARC->n_atom,1/3.0); +>>>>>>> variant B + // *avgdis = pow((pSPARC->range_x * pSPARC->range_y * pSPARC->range_z)/pSPARC->n_atom,1/3.0); +======= end int atm2, ityp2, cc2, pairIndex; cc = 0; for(ityp = 0; ityp < pSPARC->Ntypes; ityp++){ @@ -2284,7 +2490,13 @@ void PrintMD(SPARC_OBJ *pSPARC, int Flag, int print_restart_typ) { fprintf(mdout,":CELL: %.15g %.15g %.15g\n",pSPARC->range_x,pSPARC->range_y,pSPARC->range_z); //(no variable for position of barostat variable) else fprintf(mdout,":LATVEC_SCALE: %.15g %.15g %.15g\n",pSPARC->range_x/pSPARC->initialLatVecLength[0],pSPARC->range_y/pSPARC->initialLatVecLength[1],pSPARC->range_z/pSPARC->initialLatVecLength[2]); +<<<<<<< variant A fprintf(mdout,":TTHRMI(K): %.15g\n", pSPARC->thermos_Ti); +>>>>>>> variant B + if (pSPARC->CyclixFlag) + fprintf(mdout,":TWIST_ANGLE: %.15g\n", pSPARC->twist); + fprintf(mdout,":TTHRMI(K): %.15g\n", pSPARC->thermos_Ti); +======= end fprintf(mdout,":TARGET_PRESSURE: %.15g GPa\n",pSPARC->prtarget * 29421.02648438959); fprintf(mdout,":NPT_NP_ini_Hamiltonian: %.15g\n", pSPARC->init_Hamil_NPT_NP); } @@ -2335,7 +2547,11 @@ void RestartMD(SPARC_OBJ *pSPARC) { l_buff = (2 + 1) * sizeof(int) + (6 * pSPARC->n_atom + (5 + 3*pSPARC->NPT_NHnnos + 9)) * sizeof(double); } else if (strcmpi(pSPARC->MDMeth,"NPT_NP") == 0){ +<<<<<<< variant A l_buff = 2 * sizeof(int) + (6 * pSPARC->n_atom + (5 + 14)) * sizeof(double); +>>>>>>> variant B + l_buff = 2 * sizeof(int) + (6 * pSPARC->n_atom + (5 + 15)) * sizeof(double); +======= end } else { l_buff = 2 * sizeof(int) + (6 * pSPARC->n_atom + 5) * sizeof(double); @@ -2485,6 +2701,11 @@ void RestartMD(SPARC_OBJ *pSPARC) { } else if (strcmpi(str,":TTHRMI(K):") == 0) fscanf(rst_fp,"%lf", &pSPARC->thermos_Ti); +<<<<<<< variant A +>>>>>>> variant B + else if (strcmpi(str,":TWIST_ANGLE:") == 0) + fscanf(rst_fp,"%lf", &pSPARC->twist); +======= end else if (strcmpi(str,":TARGET_PRESSURE:") == 0) fscanf(rst_fp,"%lf", &pSPARC->prtarget); else if (strcmpi(str,":NPT_NP_ini_Hamiltonian:") == 0) @@ -2507,7 +2728,11 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Pack(&pSPARC->xi_nose, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->thermos_Ti, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); } +<<<<<<< variant A else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ +>>>>>>> variant B + else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ +======= end MPI_Pack(&pSPARC->NPT_NHnnos, 1, MPI_INT, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(pSPARC->NPT_NHqmass, pSPARC->NPT_NHnnos, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(pSPARC->vlogs, pSPARC->NPT_NHnnos, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); @@ -2535,6 +2760,10 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Pack(&pSPARC->range_y, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->range_z, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->thermos_Ti, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); +<<<<<<< variant A +>>>>>>> variant B + MPI_Pack(&pSPARC->twist, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); +======= end MPI_Pack(&pSPARC->prtarget, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->init_Hamil_NPT_NP, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); } @@ -2566,7 +2795,11 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Unpack(buff, l_buff, &position, &pSPARC->xi_nose, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->thermos_Ti, 1, MPI_DOUBLE, MPI_COMM_WORLD); } +<<<<<<< variant A else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ +>>>>>>> variant B + else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ +======= end MPI_Unpack(buff, l_buff, &position, &pSPARC->NPT_NHnnos, 1, MPI_INT, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, pSPARC->NPT_NHqmass, pSPARC->NPT_NHnnos, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, pSPARC->vlogs, pSPARC->NPT_NHnnos, MPI_DOUBLE, MPI_COMM_WORLD); @@ -2595,6 +2828,10 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Unpack(buff, l_buff, &position, &pSPARC->range_y, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->range_z, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->thermos_Ti, 1, MPI_DOUBLE, MPI_COMM_WORLD); +<<<<<<< variant A +>>>>>>> variant B + MPI_Unpack(buff, l_buff, &position, &pSPARC->twist, 1, MPI_DOUBLE, MPI_COMM_WORLD); +======= end MPI_Unpack(buff, l_buff, &position, &pSPARC->prtarget, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->init_Hamil_NPT_NP, 1, MPI_DOUBLE, MPI_COMM_WORLD); } @@ -2619,4 +2856,7 @@ void Rename_restart(SPARC_OBJ *pSPARC) { if( access(pSPARC->restartC_Filename, F_OK ) != -1 ) rename(pSPARC->restartC_Filename, pSPARC->restartP_Filename); } +<<<<<<< variant A +>>>>>>> variant B +======= end diff --git a/src/relax.c b/src/relax.c index d1770b91..59c53157 100644 --- a/src/relax.c +++ b/src/relax.c @@ -107,7 +107,8 @@ void NLCG(SPARC_OBJ *pSPARC) { Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*atm], &pSPARC->atom_pos[3*atm+1], &pSPARC->atom_pos[3*atm+2]); } } - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); err = delta_new = 0.0; for(atmc = 0; atmc < szatm; atmc++){ r[atmc] = pSPARC->forces[atmc]; @@ -117,7 +118,8 @@ void NLCG(SPARC_OBJ *pSPARC) { err = fabs(pSPARC->forces[atmc]); // defined as supremum norm of force vector } } else { - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->d = (double *)malloc(szatm * sizeof(double)); // search direction if (pSPARC->d == NULL) { printf("\nCannot allocate memory for search direction array in NLCG!\n"); @@ -191,9 +193,11 @@ void NLCG(SPARC_OBJ *pSPARC) { // Charge extrapolation (for better rho_guess) pSPARC->Relax_fac = sigma; + //elecDensExtrapolation converts to cell coords if cell_typ != 0 and MDFlag or RelaxFlag == 1 elecDensExtrapolation(pSPARC); Check_atomlocation(pSPARC); - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; eta_prev = 0.0; for(atmc = 0; atmc < szatm; atmc++){ @@ -220,7 +224,8 @@ void NLCG(SPARC_OBJ *pSPARC) { pSPARC->Relax_fac = alpha; elecDensExtrapolation(pSPARC); Check_atomlocation(pSPARC); - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; eta_prev = eta; j++; @@ -342,14 +347,16 @@ void LBFGS(SPARC_OBJ *pSPARC) { Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*atm], &pSPARC->atom_pos[3*atm+1], &pSPARC->atom_pos[3*atm+2]); } } - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); err = 0.0; for(i = 0; i < n; i++){ if (fabs(pSPARC->forces[i]) > err) err = fabs(pSPARC->forces[i]); // defined as supremum norm of force vector } } else { - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); err = 0.0; for(i = 0; i < n; i++){ if (fabs(pSPARC->forces[i]) > err) @@ -630,7 +637,8 @@ void LBFGS(SPARC_OBJ *pSPARC) { elecDensExtrapolation(pSPARC); Check_atomlocation(pSPARC); - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; err = 0.0; for(i = 0; i < n; i++){ @@ -743,7 +751,8 @@ void FIRE(SPARC_OBJ *pSPARC) { Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*atm], &pSPARC->atom_pos[3*atm+1], &pSPARC->atom_pos[3*atm+2]); } } - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); err = 0.0; for(i = 0; i < n; i++){ if (fabs(pSPARC->forces[i]) > err) @@ -753,7 +762,8 @@ void FIRE(SPARC_OBJ *pSPARC) { acc[i] = pSPARC->forces[i]/mass; } } else { - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); err = 0.0; for(i = 0; i < n; i++){ if (fabs(pSPARC->forces[i]) > err) @@ -829,7 +839,8 @@ void FIRE(SPARC_OBJ *pSPARC) { pSPARC->Relax_fac = 1.0; elecDensExtrapolation(pSPARC); Check_atomlocation(pSPARC); - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; for(i = 0; i < n; i++){ @@ -978,7 +989,8 @@ double volrelax_constraint(SPARC_OBJ *pSPARC, double vol) if (pSPARC->RelaxFlag == 2) { // turn stress calculation on pSPARC->Calc_stress = 1; - Calculate_electronicGroundState(pSPARC); + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); pSPARC->elecgs_Count++; pSPARC->RelaxCount++; diff --git a/src/socket/conf.h b/src/socket/conf.h new file mode 100644 index 00000000..86a7b27f --- /dev/null +++ b/src/socket/conf.h @@ -0,0 +1,9 @@ +/** +* @file conf.h +* @brief Copied conf file from libsocket. Please do not modify this file. +* @authors T.Tian +**/ +#define LIBSOCKET_VERSION 2.4 +#define LIBSOCKET_LINUX 1 +#define LIBSOCKET_FREEBSD 0 +#define LIBSOCKET_SUNOS 0 diff --git a/src/socket/driver.c b/src/socket/driver.c new file mode 100644 index 00000000..95929c1c --- /dev/null +++ b/src/socket/driver.c @@ -0,0 +1,1164 @@ +/** + * @file driver.c + * @brief Driver mode for SPARC (socket communication) + * This file contains both the socket communication methods and a socket-aware main loop + * @authors T.Tian + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include "driver.h" +#include "isddft.h" +#include "libinetsocket.h" +#include "libunixsocket.h" + +#include "initialization.h" +/* #include "orbitalElecDensInit.h" */ +#include "electronicGroundState.h" +#include "relax.h" +/* #include "md.h" */ + +/* Tools to be used in socket interface. Don't need #ifdef switches */ + +int split_socket_name(const char *str, char *host, int *port, int *inet) +{ + /** + * @brief split socket name by delimiter ":" + * examples: + * "localhost:1234" --> host = "localhost", port = 1234, inet = 1 + * "unix_socket:UNIX" --> host = "unix_socket", port = 0, inet = 0 + * ":1234" --> host = "localhost", port = 1234, inet = 1 + * "unix_socket:" --> invalid + * "localhost:" --> invalid + * "unix_socket:UNIX" + * @param str: socket name + * @param host: host name + * @param port: port number + * @param inet: 1 for inet socket, 0 for unix socket + **/ + const char *delim_pos = strchr(str, ':'); + if (delim_pos == NULL || strchr(delim_pos + 1, ':') != NULL) + { + printf("Error: string must have format host:port or unix_socket:UNIX\n"); + return 1; + } + + int delim_index = delim_pos - str; + + if (delim_index > 0) + { + strncpy(host, str, delim_index); + host[delim_index] = '\0'; + } + else + { + host[0] = '\0'; // Empty string + } + /* If the length of host is zero, then use localhost instead */ + if (strlen(host) == 0) + { + strcpy(host, "localhost"); + } + + const char *port_str = &str[delim_index + 1]; + char *endptr; + long int value = strtol(port_str, &endptr, 10); + if (*endptr == '\0') + { + // The port string was successfully converted to an integer + // value == 0 may come from "localhost:0" or empty string, + // both are not allowed + if (value > 0 && value <= 65535) + { + *port = (int)value; + } + else + { + printf("Error: port number must be a positive integer no larger than 65535\n"); + return 1; + } + } + else if (strcasecmp(port_str, "UNIX") == 0) + { + // The port string is "UNIX" (case insensitive) + *port = 0; + } + else + { + // The port string is neither a valid integer nor "UNIX" + printf("Error: port must be an integer or 'UNIX'\n"); + return 1; + } + + if (*port == 0) + { + *inet = 0; + } + else + { + *inet = 1; + } + + return 0; +} + +int initialize_Socket(SPARC_OBJ *pSPARC) +{ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + if (pSPARC->SocketFlag != 1) + { + if (rank == 0) + { + printf("You're using initializing the socket with SocketFlag != 0! Exit"); + } + exit(EXIT_FAILURE); + } + // Initialize the flags compatible with SocketFlag + // TODO: maybe we should move to another location, like copy_sparc_input + pSPARC->MDFlag = 0; + pSPARC->RelaxFlag = 0; + pSPARC->PrintAtomPosFlag = 1; + pSPARC->PrintForceFlag = 1; + pSPARC->SocketSCFCount = 0; + +#ifdef DEBUG + if (rank == 0) + printf("##########################Initializing socket##########################\n"); +#endif + // Create a fd according to the socket type + int socket_fd = -1; + // For some reason the inet socket may die out earlier. For first stage testing use unix socket only + // TODO: should fix this later when server side code is done + if (rank != 0) + { + return 0; + } + + // Only init on rank 0 + if (pSPARC->socket_inet == 1) + { + // Convert port number to string + char service[L_STRING]; + sprintf(service, "%d", pSPARC->socket_port); +#ifdef DEBUG + printf("This is your service: %s\n", service); + printf("This is your service: %s\n", pSPARC->socket_host); +#endif // DEBUG + socket_fd = create_inet_stream_socket(pSPARC->socket_host, service, LIBSOCKET_IPv4, 0); + if (socket_fd == -1) + { + printf("Failed to create or communicate with the inet socket %s:%s, Exiting...\n", + pSPARC->socket_host, service); + exit(EXIT_FAILURE); + } + pSPARC->socket_fd = socket_fd; + } + else if (pSPARC->socket_inet == 0) + { + // Use unix socket + // TODO: change the actual unix socket path under /tmp/ + // TODO: unix socket must have server side up, otherwise will die immediately + socket_fd = create_unix_stream_socket(pSPARC->socket_host, 0); + if (socket_fd == -1) + { + printf("Failed to create or communicate with the unix socket %s, Exiting...\n", + pSPARC->socket_host); + exit(EXIT_FAILURE); + } + pSPARC->socket_fd = socket_fd; + } + else + { + printf("Incorrect socket type, exiting...\n"); + printf("The code is buggy! Exiting...\n"); + exit(EXIT_FAILURE); + } +#ifdef DEBUG + printf("Created socket fd %i\n", pSPARC->socket_fd); +#endif + return 0; +} + +/** + * TODO: complete + */ +int close_Socket(SPARC_OBJ *pSPARC) +{ + int socket_fd = pSPARC->socket_fd; + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + if (rank != 0) + { + return 0; + } + + // rank 0 + int close_res = close(socket_fd); + if (close_res != 0) + { + perror("Failed to close socket!"); + exit(EXIT_FAILURE); + } + else + { +#ifdef DEBUG + printf("@Driver mode: closed socker %d.\n", pSPARC->socket_fd); +#endif + } + pSPARC->socket_fd = -1; + return 0; +} + +/** + *@brief General purpose read function. Will read until the fd until buffer len is filled + **/ +int readBuffer(int fd, void *data, int len) +{ + int n, nr; + printf("Waiting to read.......\n"); + n = nr = read(fd, data, len); + printf("After reading, there are n %d\n", n); + while (nr > 0 && n < len) + { + printf("Read %d bytes", nr); + nr = read(fd, data + n, len - n); + n += nr; + } + if (n == 0) + { + perror("Error reading from socket: server has quit or connection broke"); + return -1; + } + return 0; +} + +/** + * TODO: complete + */ +// int writeBuffer(SPARC_OBJ *pSPARC, const char *data, int len){ +// return 0; +// } + +/** + * @brief Write functions + All writeBuffer_* functions are meant to be used on rank == 0! + **/ +int writeBuffer_double(SPARC_OBJ *pSPARC, const double data_d) +{ + int sockfd = pSPARC->socket_fd; + int ret = write(sockfd, &data_d, sizeof(double)); + if (ret < 0) + fprintf(stderr, "Failed to write to socket fd %d, exiting...\n", sockfd); + return ret; +} + +int writeBuffer_int(SPARC_OBJ *pSPARC, const int data_i) +{ + int sockfd = pSPARC->socket_fd; + int ret = write(sockfd, &data_i, sizeof(int)); + if (ret < 0) + fprintf(stderr, "Failed to write to socket fd %d, exiting...\n", sockfd); + return ret; +} + +int writeBuffer_char(SPARC_OBJ *pSPARC, const char c) +{ + int sockfd = pSPARC->socket_fd; + int ret = write(sockfd, &c, sizeof(char)); + if (ret < 0) + fprintf(stderr, "Failed to write to socket fd %d, exiting...\n", sockfd); + return ret; +} + +int writeBuffer_double_vec(SPARC_OBJ *pSPARC, const double *data_dv, int len) +{ + int sockfd = pSPARC->socket_fd; + int ret = write(sockfd, data_dv, sizeof(double) * len); + if (ret < 0) + fprintf(stderr, "Failed to write to socket fd %d, exiting...\n", sockfd); + return ret; +} + +int writeBuffer_int_vec(SPARC_OBJ *pSPARC, const int *data_iv, int len) +{ + int sockfd = pSPARC->socket_fd; + int ret = write(sockfd, data_iv, sizeof(int) * len); + if (ret < 0) + fprintf(stderr, "Failed to write to socket fd %d, exiting...\n", sockfd); + return ret; +} + +int writeBuffer_string(SPARC_OBJ *pSPARC, const char *str, int len) +{ + int sockfd = pSPARC->socket_fd; + int ret = write(sockfd, str, sizeof(char) * len); + if (ret < 0) + fprintf(stderr, "Failed to write to socket fd %d, exiting...\n", sockfd); + return ret; +} + +/** + * @brief Read functions. All readBuffer_* functions are meant to be used + at rank == 0 + * **/ + +int readBuffer_double(SPARC_OBJ *pSPARC, double *data_d) +{ + int sockfd = pSPARC->socket_fd; + int ret = readBuffer(sockfd, data_d, sizeof(double)); + return ret; +} + +int readBuffer_int(SPARC_OBJ *pSPARC, int *data_i) +{ + int sockfd = pSPARC->socket_fd; + int ret = readBuffer(sockfd, data_i, sizeof(int)); + return ret; +} + +int readBuffer_char(SPARC_OBJ *pSPARC, char *c) +{ + int sockfd = pSPARC->socket_fd; + int ret = readBuffer(sockfd, c, sizeof(char)); + return ret; +} + +int readBuffer_double_vec(SPARC_OBJ *pSPARC, double *data_dv, int len) +{ + int sockfd = pSPARC->socket_fd; + int ret = readBuffer(sockfd, data_dv, sizeof(double) * len); + return ret; +} + +int readBuffer_int_vec(SPARC_OBJ *pSPARC, int *data_iv, int len) +{ + int sockfd = pSPARC->socket_fd; + int ret = readBuffer(sockfd, data_iv, sizeof(int) * len); + return ret; +} + +int readBuffer_string(SPARC_OBJ *pSPARC, char *str, int len) +{ + int sockfd = pSPARC->socket_fd; + int ret = readBuffer(sockfd, str, sizeof(char) * len); + return ret; +} + +/**@brief Function to reinitialize the mesh for SPARC object after changing the lattice + pSPARC->LatVec should already be assigned + **/ +void reinit_mesh(SPARC_OBJ *pSPARC){ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + double range_x, range_y, range_z; + /* double old_range_x, old_range_y, old_range_z; */ + + // TODO: should we disable latvec_scale in socket mode? + range_x = sqrt(pSPARC->LatVec[0] * pSPARC->LatVec[0] + pSPARC->LatVec[1] * pSPARC->LatVec[1] + pSPARC->LatVec[2] * pSPARC->LatVec[2]); + range_y = sqrt(pSPARC->LatVec[3] * pSPARC->LatVec[3] + pSPARC->LatVec[4] * pSPARC->LatVec[4] + pSPARC->LatVec[5] * pSPARC->LatVec[5]); + range_z = sqrt(pSPARC->LatVec[6] * pSPARC->LatVec[6] + pSPARC->LatVec[7] * pSPARC->LatVec[7] + pSPARC->LatVec[8] * pSPARC->LatVec[8]); + /* old_range_x = pSPARC->range_x; */ + pSPARC->range_x = range_x; + /* old_range_y = pSPARC->range_y; */ + pSPARC->range_y = range_y; + /* old_range_z = pSPARC->range_z; */ + pSPARC->range_z = range_z; + + // Update the lattice matrix information like Jacbdet, LatUVec, etc + Cart2nonCart_transformMat(pSPARC); + + + // Update finite difference delta + pSPARC->delta_x = pSPARC->range_x/(pSPARC->numIntervals_x); + pSPARC->delta_y = pSPARC->range_y/(pSPARC->numIntervals_y); + pSPARC->delta_z = pSPARC->range_z/(pSPARC->numIntervals_z); + pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z * pSPARC->Jacbdet; + + int FDn = pSPARC->order / 2; + int p, i; + // 1st derivative weights including mesh + double dx_inv, dy_inv, dz_inv; + dx_inv = 1.0 / pSPARC->delta_x; + dy_inv = 1.0 / pSPARC->delta_y; + dz_inv = 1.0 / pSPARC->delta_z; + for (p = 1; p < FDn + 1; p++) { + pSPARC->D1_stencil_coeffs_x[p] = pSPARC->FDweights_D1[p] * dx_inv; + pSPARC->D1_stencil_coeffs_y[p] = pSPARC->FDweights_D1[p] * dy_inv; + pSPARC->D1_stencil_coeffs_z[p] = pSPARC->FDweights_D1[p] * dz_inv; + } + + + // 2nd derivative weights including mesh + double dx2_inv, dy2_inv, dz2_inv; + dx2_inv = 1.0 / (pSPARC->delta_x * pSPARC->delta_x); + dy2_inv = 1.0 / (pSPARC->delta_y * pSPARC->delta_y); + dz2_inv = 1.0 / (pSPARC->delta_z * pSPARC->delta_z); + + // Stencil coefficients for mixed derivatives + if (pSPARC->cell_typ == 0) { + for (p = 0; p < FDn + 1; p++) { + pSPARC->D2_stencil_coeffs_x[p] = pSPARC->FDweights_D2[p] * dx2_inv; + pSPARC->D2_stencil_coeffs_y[p] = pSPARC->FDweights_D2[p] * dy2_inv; + pSPARC->D2_stencil_coeffs_z[p] = pSPARC->FDweights_D2[p] * dz2_inv; + } + } else if (pSPARC->cell_typ > 10 && pSPARC->cell_typ < 20) { + for (p = 0; p < FDn + 1; p++) { + pSPARC->D2_stencil_coeffs_x[p] = pSPARC->lapcT[0] * pSPARC->FDweights_D2[p] * dx2_inv; + pSPARC->D2_stencil_coeffs_y[p] = pSPARC->lapcT[4] * pSPARC->FDweights_D2[p] * dy2_inv; + pSPARC->D2_stencil_coeffs_z[p] = pSPARC->lapcT[8] * pSPARC->FDweights_D2[p] * dz2_inv; + pSPARC->D2_stencil_coeffs_xy[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dx_inv; // 2*T_12 d/dx(df/dy) + pSPARC->D2_stencil_coeffs_xz[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dx_inv; // 2*T_13 d/dx(df/dz) + pSPARC->D2_stencil_coeffs_yz[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // 2*T_23 d/dy(df/dz) + pSPARC->D1_stencil_coeffs_xy[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dy_inv; // d/dx(2*T_12 df/dy) used in d/dx(2*T_12 df/dy + 2*T_13 df/dz) + pSPARC->D1_stencil_coeffs_yx[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dx_inv; // d/dy(2*T_12 df/dx) used in d/dy(2*T_12 df/dx + 2*T_23 df/dz) + pSPARC->D1_stencil_coeffs_xz[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dz_inv; // d/dx(2*T_13 df/dz) used in d/dx(2*T_12 df/dy + 2*T_13 df/dz) + pSPARC->D1_stencil_coeffs_zx[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dx_inv; // d/dz(2*T_13 df/dx) used in d/dz(2*T_13 df/dz + 2*T_23 df/dy) + pSPARC->D1_stencil_coeffs_yz[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dz_inv; // d/dy(2*T_23 df/dz) used in d/dy(2*T_12 df/dx + 2*T_23 df/dz) + pSPARC->D1_stencil_coeffs_zy[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // d/dz(2*T_23 df/dy) used in d/dz(2*T_12 df/dx + 2*T_23 df/dy) + } + } + + // maximum eigenvalue of -0.5 * Lap (only with periodic boundary conditions) + if(pSPARC->cell_typ == 0) { + pSPARC->MaxEigVal_mhalfLap = pSPARC->D2_stencil_coeffs_x[0] + + pSPARC->D2_stencil_coeffs_y[0] + + pSPARC->D2_stencil_coeffs_z[0]; + double scal_x, scal_y, scal_z; + scal_x = (pSPARC->Nx - pSPARC->Nx % 2) / (double) pSPARC->Nx; + scal_y = (pSPARC->Ny - pSPARC->Ny % 2) / (double) pSPARC->Ny; + scal_z = (pSPARC->Nz - pSPARC->Nz % 2) / (double) pSPARC->Nz; + for (int p = 1; p < FDn + 1; p++) { + pSPARC->MaxEigVal_mhalfLap += 2.0 * (pSPARC->D2_stencil_coeffs_x[p] * cos(M_PI*p*scal_x) + + pSPARC->D2_stencil_coeffs_y[p] * cos(M_PI*p*scal_y) + + pSPARC->D2_stencil_coeffs_z[p] * cos(M_PI*p*scal_z)); + } + pSPARC->MaxEigVal_mhalfLap *= -0.5; + } + + double h_eff = 0.0; + if (fabs(pSPARC->delta_x - pSPARC->delta_y) < 1E-12 && + fabs(pSPARC->delta_y - pSPARC->delta_z) < 1E-12) { + h_eff = pSPARC->delta_x; + } else { + // find effective mesh s.t. it has same spectral width + h_eff = sqrt(3.0 / (dx2_inv + dy2_inv + dz2_inv)); + } + + // find Chebyshev polynomial degree based on max eigenvalue (spectral width) + if (pSPARC->ChebDegree < 0) { + pSPARC->ChebDegree = Mesh2ChebDegree(h_eff); + } else { + } + + // default Kerker tolerance + if (pSPARC->TOL_PRECOND < 0.0) { // kerker tol not provided by user + pSPARC->TOL_PRECOND = (h_eff * h_eff) * 1e-3; + } + + + Calculate_kpoints(pSPARC); + if (pSPARC->Nkpts >= 1 && pSPARC->kptcomm_index != -1) { + Calculate_local_kpoints(pSPARC); + } + Calculate_PseudochargeCutoff(pSPARC); + + + if (rank == 0){ + write_output_reinit(pSPARC); + } +} + +/** @brief Function to map the atomic coordinates to the domain within LatVec. + * will throw error if the positions are outside domain + **/ +void map_atom_coord(SPARC_OBJ *pSPARC){ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + int i = 0; + // map positions back into the domain if necessary, the function is taken from initialization.c + if (pSPARC->BCx == 0) { + double x; + for (i = 0; i < pSPARC->n_atom; i++) { + x = *(pSPARC->atom_pos+3*i); + if (x < 0 || x > pSPARC->range_x) + { + x = fmod(x,pSPARC->range_x); + *(pSPARC->atom_pos+3*i) = x + (x<0.0)*pSPARC->range_x; + } + } + } else if (rank == 0) { + // for Dirichlet BC, exit if atoms goes out of the domain + double x; + for (i = 0; i < pSPARC->n_atom; i++) { + x = *(pSPARC->atom_pos+3*i); + if (x < pSPARC->xin || x > pSPARC->range_x + pSPARC->xin) + { + printf("\nERROR: position of atom # %d is out of the domain with Dirichlet BC!\n",i+1); + exit(EXIT_FAILURE); + } + } + } + + if (pSPARC->BCy == 0) { + double y; + for (i = 0; i < pSPARC->n_atom; i++) { + y = *(pSPARC->atom_pos+1+3*i); + if (y < 0 || y > pSPARC->range_y) + { + y = fmod(y,pSPARC->range_y); + *(pSPARC->atom_pos+1+3*i) = y + (y<0.0)*pSPARC->range_y; + } + } + } else if (rank == 0) { + // for Dirichlet BC, exit if atoms goes out of the domain + double y; + for (i = 0; i < pSPARC->n_atom; i++) { + y = *(pSPARC->atom_pos+1+3*i); + if (y < 0 || y > pSPARC->range_y) + { + printf("\nERROR: position of atom # %d is out of the domain with Dirichlet BC!\n",i+1); + exit(EXIT_FAILURE); + } + } + } + + if (pSPARC->BCz == 0) { + double z; + for (i = 0; i < pSPARC->n_atom; i++) { + z = *(pSPARC->atom_pos+2+3*i); + if (z < 0 || z > pSPARC->range_z) + { + z = fmod(z,pSPARC->range_z); + *(pSPARC->atom_pos+2+3*i) = z + (z<0.0)*pSPARC->range_z; + } + } + } else if (rank == 0){ + // for Dirichlet BC, exit if atoms goes out of the domain + double z; + for (i = 0; i < pSPARC->n_atom; i++) { + z = *(pSPARC->atom_pos+2+3*i); + if (z < 0 || z > pSPARC->range_z) + { + printf("\nERROR: position of atom # %d is out of the domain with Dirichlet BC!\n",i+1); + exit(EXIT_FAILURE); + } + } + } +} + +/** + * @brief Function that reassign atomic positions / lattice vectors etc to SPARC_OBJ + * This function will be used at each initialization / single point loop + * This function is MPI safe + * TODO: check if need strict atoms size check + * **/ + + + +void reassign_atoms_info(SPARC_OBJ *pSPARC, int natoms, double *atom_pos, double *lattice, double *reci_lattice) +{ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + // For now only support natoms == n_atom + if (natoms != pSPARC->n_atom) + { + if (rank == 0) + printf("SPARC's socket interface currently does not support changing number of atoms yet! Exit...\n"); + exit(EXIT_FAILURE); + } + + // TODO: fix the cell type check + + + pSPARC->atom_pos = (double *)malloc(sizeof(double) * 3 * natoms); + memcpy(pSPARC->atom_pos, atom_pos, sizeof(double) * 3 * natoms); + memcpy(pSPARC->LatVec, lattice, sizeof(double) * 9); + map_atom_coord(pSPARC); + reinit_mesh(pSPARC); +} + +/** + * @brief Read message from sparc and return a status code + * Use status defined in socket.h. + This is a function for ALL RANKS + **/ +int read_socket_header(SPARC_OBJ *pSPARC, int *status) +{ + int rank, size, rank0_finish; + double sleep_interval = 0.1; // Checking every 0.5 seconds + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + rank0_finish = 0; // Helper for other ranks to check + int ret = 0; + char header[IPI_HEADERLEN]; + int sockfd = pSPARC->socket_fd; + MPI_Request request; + + if (rank == 0) + { + ret = readBuffer_string(pSPARC, header, IPI_HEADERLEN); + rank0_finish = 1; +#ifdef DEBUG + printf("@Driver mode: get raw header %s\n", header); +#endif + // Status code should broadcast to + if (ret < 0) + { + fprintf(stderr, "Failed to read from socket fd %d, exiting...\n", sockfd); + } + // Trim the string and compare with upper case letters + else if (strncasecmp(header, "STATUS", strlen("STATUS")) == 0) + { + *status = IPI_MSG_STATUS; + } + else if (strncasecmp(header, "POSDATA", strlen("POSDATA")) == 0) + { + *status = IPI_MSG_POSDATA; + } + else if (strncasecmp(header, "GETFORCE", strlen("GETFORCE")) == 0) + { + *status = IPI_MSG_GETFORCE; + } + /* Extra keywords for extended SPARC protocol */ + else if (strncasecmp(header, "SETPARAM", strlen("SETPARAM")) == 0) + { + *status = SPARC_MSG_SETPARAM; + } + else + { + *status = IPI_MSG_OTHER; + } + /* Notify other ranks */ + for (int i = 1; i < size; i++){ + MPI_Isend(&rank0_finish, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &request); + } + } + else { + MPI_Irecv(&rank0_finish, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &request); + int ready = 0; + while (rank0_finish == 0){ + MPI_Test(&request, &ready, MPI_STATUS_IGNORE); + if (ready == 0) + sleep(sleep_interval); + else if (rank0_finish == 0){ + printf("Received rank0_finish %d at rank %d, this is a bug.\n", rank0_finish, rank); + exit(1); + } + } + } + MPI_Barrier(MPI_COMM_WORLD); + MPI_Bcast(&ret, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(status, 1, MPI_INT, 0, MPI_COMM_WORLD); + return ret; +} + +/** + * @brief Get positions from socket and reassign to SPARC_OBJ + * MPI broadcast is done in this function + This is a function for ALL RANKS + * + **/ +int read_atoms_position_fom_socket(SPARC_OBJ *pSPARC, int init) +{ + int rank; + int ret = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + // TODO: check if status is in IPI_MSG_POSDATA + int sockfd = pSPARC->socket_fd; + int natoms; + double ipi_cell[9], ipi_inv_cell[9]; + // The read process will need to get the forces, although we don't pass them back to SPARC + double *ipi_atom_pos, *ipi_forces; + // TODO: do read at rank0 + if (rank == 0){ + ret |= readBuffer_double_vec(pSPARC, ipi_cell, 9); + ret |= readBuffer_double_vec(pSPARC, ipi_inv_cell, 9); + ret |= readBuffer_int(pSPARC, &natoms); + } + MPI_Bcast(&ret, 1, MPI_INT, 0, MPI_COMM_WORLD); + if (ret < 0){ + perror("Error receiving Cell and Natoms information from socket! Will exit"); + return ret; + } + // Bcast info to each rank + MPI_Bcast(&natoms, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(ipi_cell, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(ipi_inv_cell, 9, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + // on each rank + ipi_atom_pos = (double *)malloc(sizeof(double) * 3 * natoms); + ipi_forces = (double *)malloc(sizeof(double) * 3 * natoms); + + // Read real positions, now rank 0 + if (rank == 0) + { + ret |= readBuffer_double_vec(pSPARC, ipi_atom_pos, 3 * natoms); + } + MPI_Bcast(&ret, 1, MPI_INT, 0, MPI_COMM_WORLD); + if (ret < 0){ + perror("Error receiving Pos information from socket! Will exit"); + return ret; + } + + MPI_Bcast(ipi_atom_pos, 3 * natoms, MPI_DOUBLE, 0, MPI_COMM_WORLD); + // TODO: make sure the values are broadcast + // TODO: wrap cell for the positions +#ifdef DEBUG + if (rank == 0) + { + printf("Starting socket communication\n"); + printf("Received from socket the following position data:\n"); + printf("natoms: %d\n", natoms); + printf("cell: \t%f %f %f \n\t%f %f %f \n\t%f %f %f\n", ipi_cell[0], ipi_cell[1], ipi_cell[2], ipi_cell[3], ipi_cell[4], ipi_cell[5], ipi_cell[6], ipi_cell[7], ipi_cell[8]); + printf("inverse cell \t%f %f %f \n\t%f %f %f \n\t%f %f %f\n", ipi_inv_cell[0], ipi_inv_cell[1], ipi_inv_cell[2], ipi_inv_cell[3], ipi_inv_cell[4], ipi_inv_cell[5], ipi_inv_cell[6], ipi_inv_cell[7], ipi_inv_cell[8]); + } +#endif // DEBUG + + // This should work on each rank + reassign_atoms_info(pSPARC, natoms, ipi_atom_pos, ipi_cell, ipi_inv_cell); + + free(ipi_atom_pos); + free(ipi_forces); + return ret; +} + +/** + * @brief Convert stress 6 vector to 9 matrix virial + * SPARC's stress information (in PBC) is [xx, xy, xz, yy, yz, zz] in Ha/Bohr^3, different from the common voigt stress + * Virial will be [xx, xy, xz, xy, yy, yz, xz, yz, zz] * -volume in Ha + * Cell volume is given by: pSPARC->Jacbdet * pSPARC->range_x * pSPARC->range_y * pSPARC->range_z + **/ +void stress_to_virial(SPARC_OBJ *pSPARC, double *virial) +{ + // TODO: may need to be a bit looseen + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + double volCell = pSPARC->Jacbdet * pSPARC->range_x * pSPARC->range_y * pSPARC->range_z; + + // This is the general form for the converted virial notation of the stress. For PBC (pSPARC->BC == 2) it's Ha/Bohr^3 + // In other BCs, we should do post processing + double virial_calc[9] = { + -pSPARC->stress[0] * volCell, -pSPARC->stress[1] * volCell, -pSPARC->stress[2] * volCell, + -pSPARC->stress[1] * volCell, -pSPARC->stress[3] * volCell, -pSPARC->stress[4] * volCell, + -pSPARC->stress[2] * volCell, -pSPARC->stress[4] * volCell, -pSPARC->stress[5] * volCell + }; + + // Stress not calculated, return zero-tensor + if (pSPARC->Calc_stress == 0){ + if (rank == 0) + printf("Stress is not calculated in this system. Will return zero-stress tensor to the socket server.\n"); + for (int i; i < 9; i++){ + virial_calc[i]= 0.0; + } + } + // 0D case. make sure all components are 0 + else if (pSPARC->BC == 1){ + if (rank == 0) + printf("You're simulating an isolated system, no stress will be calculated. Zero stress tensor will be send to the socket server\n"); + for (int i; i < 9; i++){ + virial_calc[i]= 0.0; + } + } + // 2D stress cases + else if (pSPARC->BC == 3){ + if (rank == 0) + printf("The simulation has 2D stress tensor in Ha/Bohr**2. The equivalent value in periodic system will be returned!\n"); + // Periodic in xy- + if (pSPARC->BCx == 0 && pSPARC->BCy == 0){ + for (int i; i < 9; i++){ + virial_calc[i] /= pSPARC->range_z; + }} + // Periodic in xz- + else if (pSPARC->BCx == 0 && pSPARC->BCz == 0){ + for (int i; i < 9; i++){ + virial_calc[i] /= pSPARC->range_y; + }} + // Periodic in yz- + else if (pSPARC->BCy == 0 && pSPARC->BCz == 0){ + for (int i; i < 9; i++){ + virial_calc[i] /= pSPARC->range_x; + } + } + } + // 1D stress cases + else if (pSPARC->BC == 4){ + if (rank == 0) + printf("The simulation has 1D stress tensor in Ha/Bohr. The equivalent value in periodic system will be returned!\n"); + // Periodic in x- + if (pSPARC->BCx == 0){ + for (int i; i < 9; i++){ + virial_calc[i] /= (pSPARC->range_y * pSPARC->range_z); + } + } + // Periodic in y- + else if (pSPARC->BCy == 0){ + for (int i; i < 9; i++){ + virial_calc[i] /= (pSPARC->range_x * pSPARC->range_z); + } + } + else if (pSPARC->BCx == 0){ + for (int i; i < 9; i++){ + virial_calc[i] /= (pSPARC->range_x * pSPARC->range_y); + } + } + } + /* Cyclic or helix systems, copy from stress.c */ + else if (pSPARC->BC >= 5 && pSPARC->BC <= 7){ + if (rank == 0) + printf("The simulation has 1D stress tensor in Ha/Bohr (cyclic or helix boundaries). The equivalent value in periodic system will be returned!\n"); + for (int i; i < 9; i++){ + virial_calc[i] /= (M_PI * ( pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0) ) ) ; + } + } +#ifdef DEBUG + if (rank == 0) + { + printf("SPARC's electronic stress information is (Ha/Bohr^3): %f %f %f %f %f %f\n", pSPARC->stress[0], pSPARC->stress[1], pSPARC->stress[2], pSPARC->stress[3], pSPARC->stress[4], pSPARC->stress[5]); + printf("Virial matrix is (Ha): \t%f %f %f \n\t%f %f %f \n\t%f %f %f\n", virial_calc[0], virial_calc[1], virial_calc[2], virial_calc[3], virial_calc[4], virial_calc[5], virial_calc[6], virial_calc[7], virial_calc[8]); + } +#endif // DEBUG + memcpy(virial, virial_calc, sizeof(double) * 9); +} + +int write_forces_to_socket(SPARC_OBJ *pSPARC) +{ + // TODO: check if status is in IPI_MSG_POSDATA + int rank; + int ret = 1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + int sockfd = pSPARC->socket_fd; + int natoms = pSPARC->n_atom; + double ipi_potential; + double ipi_virial[9]; + // The read process will need to get the forces, although we don't pass them back to SPARC + double ipi_atom_pos[3 * natoms], ipi_forces[3 * natoms]; + // Check if we need to substract entropy + ipi_potential = pSPARC->Etot; + memcpy(ipi_forces, pSPARC->forces, sizeof(double) * 3 * natoms); + stress_to_virial(pSPARC, ipi_virial); + /* If any return value is -1, it will make final ret 0 */ + ret &= (write_message_to_socket(pSPARC, "FORCEREADY") + 1); + if (rank == 0){ + ret &= (writeBuffer_double(pSPARC, ipi_potential) + 1); + ret &= (writeBuffer_int(pSPARC, natoms) + 1); + ret &= (writeBuffer_double_vec(pSPARC, ipi_forces, 3 * natoms) + 1); + ret &= (writeBuffer_double_vec(pSPARC, ipi_virial, 9) + 1); + // No more message to send + // SPARC serialization should be handled by the SPARC protocol in Python API + ret &= (writeBuffer_int(pSPARC, 0) + 1); + } + ret -= 1; + if ((ret < 0) && (rank == 0)) + perror("Error writing potential / forces / stress to socket server! Exit...\n"); + MPI_Bcast(&ret, 1, MPI_INT, 0, MPI_COMM_WORLD); + return ret; +} + +/* Wrapper for writer message. This is a function for all ranks */ +int write_message_to_socket(SPARC_OBJ *pSPARC, char *message) +{ + int rank; + int ret = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + int sockfd = pSPARC->socket_fd; + char clean_msg[IPI_HEADERLEN + 1]; + // memset(clean_msg, ' ', IPI_HEADERLEN + 1); + if (strlen(message) > IPI_HEADERLEN) + { + fprintf(stderr, "Message header %s is too long to be sent through socket! Please contact the developers.\n", message); + exit(EXIT_FAILURE); + } + + strncpy(clean_msg, message, strlen(message)); + memset(clean_msg + strlen(message), ' ', IPI_HEADERLEN - strlen(message)); + clean_msg[IPI_HEADERLEN] = '\0'; + if (rank == 0){ + ret = writeBuffer_string(pSPARC, clean_msg, IPI_HEADERLEN); +#ifdef DEBUG + printf("@Driver mode: Sending message to socket: %s###\n", clean_msg); +#endif // DEBUG + } + MPI_Bcast(&ret, 1, MPI_INT, 0, MPI_COMM_WORLD); + return ret; +} + +/** + * + * + * @brief Print the atomic positions and lattice to the static file + * this is an updated version from write_output_init. + * we skip the first step since it's already initialized + */ +void socket_static_print_atom_pos(SPARC_OBJ *pSPARC) +{ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + if (rank != 0) + return; + +#ifdef DEBUG + printf("SocketSCFCOUNT is %d \n", pSPARC->SocketSCFCount); +#endif + // Do not write header for socket mode at init stage + if ((pSPARC->SocketSCFCount == 0) && (pSPARC->SocketFlag == 1)) + return; + + // SocketFlag = 0 but SocketSCFCount != 0 cannot happen! + if ((pSPARC->SocketSCFCount != 0) && (pSPARC->SocketFlag == 0)){ + perror("Internal error! SocketFlag = 0 but there is a non-zero SocketSCFCount, a bug may exist in the code. Please contact the developer."); + exit(-1); + } + + + if ((pSPARC->PrintAtomPosFlag == 1 || pSPARC->PrintForceFlag == 1) && pSPARC->MDFlag == 0 && pSPARC->RelaxFlag == 0) + { + FILE *static_fp = fopen(pSPARC->StaticFilename, "a"); + if (static_fp == NULL) + { + printf("\nCannot open file \"%s\"\n", pSPARC->StaticFilename); + exit(EXIT_FAILURE); + } + + // print atoms in either old static fashion or socket fashion + if (pSPARC->PrintAtomPosFlag == 1) + { + if (pSPARC->SocketFlag == 1){ + fprintf(static_fp, "***************************************************************************\n"); + fprintf(static_fp, " Atom positions (socket step %d) \n", pSPARC->SocketSCFCount); + fprintf(static_fp, "***************************************************************************\n"); + } + else{ + fprintf(static_fp, "***************************************************************************\n"); + fprintf(static_fp," Atom positions \n"); + fprintf(static_fp, "***************************************************************************\n"); + } + + int count = 0; + for (int i = 0; i < pSPARC->Ntypes; i++) + { + fprintf(static_fp, "Fractional coordinates of %s:\n", &pSPARC->atomType[L_ATMTYPE * i]); + for (int j = 0; j < pSPARC->nAtomv[i]; j++) + { + fprintf(static_fp, "%18.10f %18.10f %18.10f\n", + pSPARC->atom_pos[3 * count] / pSPARC->range_x, + pSPARC->atom_pos[3 * count + 1] / pSPARC->range_y, + pSPARC->atom_pos[3 * count + 2] / pSPARC->range_z); + count++; + } + } + + // Step 2: print the LATVEC information (required for each) + if (pSPARC->SocketFlag == 1){ + fprintf(static_fp, "Lattice (Bohr):\n"); + fprintf(static_fp, "%18.10E %18.10E %18.10E \n", pSPARC->LatVec[0], pSPARC->LatVec[1], pSPARC->LatVec[2]); + fprintf(static_fp, "%18.10E %18.10E %18.10E \n", pSPARC->LatVec[3], pSPARC->LatVec[4], pSPARC->LatVec[5]); + fprintf(static_fp, "%18.10E %18.10E %18.10E \n", pSPARC->LatVec[6], pSPARC->LatVec[7], pSPARC->LatVec[8]); + } + } + + fclose(static_fp); + } +} + +/* Print error information for premature socket closure to output file */ +void print_socket_error_info(SPARC_OBJ *pSPARC) +{ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + if (rank != 0) + return; + + FILE *output_fp = fopen(pSPARC->OutFilename,"w"); + if (output_fp == NULL) { + printf(stderr, "\nCannot open file \"%s\"\n",pSPARC->OutFilename); + exit(EXIT_FAILURE); + } + + fprintf(output_fp, "***************************************************************************\n"); + fprintf(output_fp," !!!WARNING!!! \n"); + fprintf(output_fp, "* Communication between SPARC and the socket server *\n"); + fprintf(output_fp, "* was dropped prematurely. Some calculation results *\n"); + fprintf(output_fp, "* may be incomplete. Please check your results. *\n"); + fprintf(output_fp, "***************************************************************************\n"); + fclose(output_fp); + +} + +/** + Below are functions used for extended SPARC protocol controls + Most of them are still in proof-of-concept stage, to showcase + that message passing is doable + **/ + +/**@brief Test function to set parameter in SPARC. For now, only demonstrate read and write are possible + **/ +int read_setparam(SPARC_OBJ *pSPARC) +{ + int rank; + int ret = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + // TODO: check if status is in IPI_MSG_POSDATA + int sockfd = pSPARC->socket_fd; + int msglen = 0; + char *msg; + if (rank == 0) + { + printf("Receiving param name \n"); + readBuffer_int(pSPARC, &msglen); + printf("Will receive a msg with length %d\n", msglen); + msg = (char *)malloc(sizeof(char) * (msglen + 1)); + readBuffer_string(pSPARC, msg, msglen); + msg[msglen] = '\0'; + printf("Received param %s\n", msg); + free(msg); + + printf("Receiving param value \n"); + readBuffer_int(pSPARC, &msglen); + printf("Will receive a msg with length %d\n", msglen); + msg = (char *)malloc(sizeof(char) * (msglen + 1)); + readBuffer_string(pSPARC, msg, msglen); + msg[msglen] = '\0'; + printf("Received value %s\n", msg); + free(msg); + } + // TODO: when it's a real function, make return values correct + return 0; +} + +/** End SPARC protocol functions **/ + +/** + * + * + * @brief Main function to implement socket communication. It should be independent of other main_Functions + */ +void main_Socket(SPARC_OBJ *pSPARC) +{ + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); +#ifdef DEBUG + if (rank == 0) + printf("Starting socket communication\n"); +#endif // DEBUG + int sfd = pSPARC->socket_fd; + int status, init, hasdata, retcode, print_socket_err; + status = -1; // -1: not initialized + init = 1; + hasdata = 0; + print_socket_err = 0; +#ifdef DEBUG + if (rank == 0) + printf("Status: socket_max_niter %d\n", pSPARC->socket_max_niter); +#endif // DEBUG + // TODO: option to specify N_MAXSTEPS (or directly taken from MD / relax?) + while (pSPARC->SocketSCFCount <= pSPARC->socket_max_niter) + { + retcode = read_socket_header(pSPARC, &status); + // We can safely close the socket in this case. + if (retcode != 0) break; + if (status == IPI_MSG_STATUS) + { + // TODO: may want to implement the NEEDINIT method + if (hasdata == 1) + { + retcode = write_message_to_socket(pSPARC, "HAVEDATA"); + } + else + { + retcode = write_message_to_socket(pSPARC, "READY"); + } + // For write functions, retcode >= 0 is success + if (retcode < 0) + { + print_socket_err = 1; + break; + } + } + else if (status == IPI_MSG_NEEDINIT) + { + if (rank == 0) + perror("NEEDINIT is not implemented yet!\n"); + exit(EXIT_FAILURE); + } + else if (status == IPI_MSG_POSDATA) + { + // Need to put the SocketSCFCount first due to print + pSPARC->SocketSCFCount++; + retcode = read_atoms_position_fom_socket(pSPARC, init); + if (retcode < 0) + { + print_socket_err = 1; + break; + } + socket_static_print_atom_pos(pSPARC); + if (init == 1) + { + init = 0; + } + // TODO: implement the electron density extrapolation method here + if (pSPARC->elecgs_Count > 0) + { + // elecDensExtrapolation(pSPARC); + // Check_atomlocation(pSPARC); + } + Calculate_Properties(pSPARC); + //Calculate_electronicGroundState(pSPARC); +#ifdef DEBUG + if (rank == 0) + { + printf("@Driver mode: single point #%d, Total energy: %f\n", pSPARC->elecgs_Count, pSPARC->Etot); + printf("@Driver mode: total energy in eV unit: %f\n", pSPARC->Etot * CONST_EH); + } +#endif // DEBUG + pSPARC->elecgs_Count++; + hasdata = 1; + } + else if (status == IPI_MSG_GETFORCE) + { + retcode = write_forces_to_socket(pSPARC); + if (retcode < 0) + { + print_socket_err = 1; + break; + } + hasdata = 0; + } + /* All extended SPARC protocol states should be listed here! */ + else if (status == SPARC_MSG_SETPARAM) + { + // The message should be received after status READY + retcode = read_setparam(pSPARC); + } + /* END of SPARC protocol settings */ + else if (status == IPI_MSG_OTHER) + { + if (rank == 0) + perror("Getting an unknown message from server, exiting...\n"); + exit(EXIT_FAILURE); + } + else if (status == IPI_MSG_EXIT) + { +#ifdef DEBUG + if (rank == 0) + printf("Server requesting SPARC to exit. Break the loop.\n"); +#endif // DEBUG + break; + } + } + if (print_socket_err != 0) + print_socket_error_info(pSPARC); + close_Socket(pSPARC); +} diff --git a/src/socket/driver.h b/src/socket/driver.h new file mode 100644 index 00000000..044e6fdf --- /dev/null +++ b/src/socket/driver.h @@ -0,0 +1,74 @@ +#ifndef DRIVER_H +#define DRIVER_H +#include "isddft.h" +#include "libinetsocket.h" +#include "libunixsocket.h" + +/** + * @brief split socket name by delimiter ":" + * examples: + * "localhost:1234" --> host = "localhost", port = 1234, inet = 1 + * "unix_socket:UNIX" --> host = "unix_socket", port = 0, inet = 0 + * ":1234" --> host = "localhost", port = 1234, inet = 1 + * "unix_socket:" --> invalid + * "localhost:" --> invalid + * "unix_socket:UNIX" + * @param str: socket name + * @param host: host name + * @param port: port number + * @param inet: 1 for inet socket, 0 for unix socket + **/ +int split_socket_name(const char *str, char *host, int *port, int *inet); + +/** + * @brief Initialize a socket file descriptor for communication with the client. + * the FD at pSPARC->socket_fd is created according to pSPARC->socket_inet + **/ +int initialize_Socket(SPARC_OBJ *pSPARC); + +/** + * @brief Destroy the socket file descriptor + **/ +int close_Socket(SPARC_OBJ *pSPARC); + +/** + * @brief Read buffer from the socket fd into string data and length of len + **/ +// int readBuffer(SPARC_OBJ *pSPARC, char *data, int len); + +/** + * @brief Write buffer of size len to the socket fd + **/ +// int writeBuffer(SPARC_OBJ *pSPARC, const char *data, int len); + + +/** + * @brief Main function with socket control + **/ +void main_Socket(SPARC_OBJ *pSPARC); + +// IPI constants + +#define IPI_HEADERLEN 12 +#define IPI_MSG_NEEDINIT 0 +#define IPI_MSG_STATUS 1 +#define IPI_MSG_POSDATA 2 +#define IPI_MSG_GETFORCE 3 +#define IPI_MSG_OTHER 4 + +/* + Below are extended MSG states for extended SPARC protocol. + New protocol keywords should be defined starting from 100 + */ +#define SPARC_MSG_SETPARAM 100 + +/* Exit status */ +#define IPI_MSG_EXIT 999 + +// Hartree to eV conversions +#define HARTREE_TO_EV 27.211386024367243 +#define HARTREE_BOHR_TO_EV_ANGSTROM 51.422067090480645 +#define HARTREE_BOHR3_TO_EV_ANGSTROM3 183.6315353072019 + + +#endif // DRIVER_H diff --git a/src/socket/libinetsocket.c b/src/socket/libinetsocket.c new file mode 100644 index 00000000..f68989b0 --- /dev/null +++ b/src/socket/libinetsocket.c @@ -0,0 +1,1106 @@ +#ifndef _GNU_SOURCE +#define _GNU_SOURCE +#endif + +#include "conf.h" +#define LIBSOCKET_VERSION 2.4 +#ifdef BD_ANDROID +#define LIBSOCKET_LINUX 0 +#else +#define LIBSOCKET_LINUX 1 +#endif + +#include +#include +#include // getaddrinfo() +#include // e.g. struct sockaddr_in on OpenBSD +#include +#include +#include +#include +#include +#include +#include +#include // read()/write() + +/** + * @file libinetsocket.c + * @brief A static copy of the libsocket library's inet socket functions. + * This library is only intended to for Linux platforms + * + * @authors T.Tian + * + * Original license information: + * This is the main file for libinetsocket. It contains all functions + * used to work with INET and INET6 sockets, both TCP and UDP. + */ +/** + * @addtogroup libinetsocket + * @{ + */ + +/* + The committers of the libsocket project, all rights reserved + (c) 2012, dermesser + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS “AS IS” AND ANY + EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY + DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +/* + * Structure of the functions defined here: + * + * + * + * + * + */ + +/* Macro definitions */ + +//# define VERBOSE // Write errors on stderr? + +#define LIBSOCKET_BACKLOG \ + 128 ///< Linux accepts a backlog value at listen() up to 128 + +// Symbolic macros +#define LIBSOCKET_TCP 1 ///< Protocol flag +#define LIBSOCKET_UDP 2 ///< Protocol flag + +#define LIBSOCKET_IPv4 3 ///< Address family flag +#define LIBSOCKET_IPv6 4 ///< Adress family flag +#define LIBSOCKET_BOTH \ + 5 ///< Adress family flag: what fits best (TCP/UDP or IPv4/6; delegate the + ///< decision to `getaddrinfo()`) + +#define LIBSOCKET_READ 1 ///< Flag for shutdown +#define LIBSOCKET_WRITE 2 ///< Flag for shutdown + +#define LIBSOCKET_NUMERIC \ + 1 ///< May be specified as flag for functions to signalize that the name + ///< resolution should not be performed. + +/** + * Writes an error to stderr without modifying errno. + */ +#define debug_write(str) \ + { \ + int verbose_errno_save = errno; \ + write(2, str, strlen(str)); \ + errno = verbose_errno_save; \ + } + +#ifdef __FreeBSD__ +#define _TRADITIONAL_RDNS +#endif + +/** + * @brief Checks return value for error. + * + * Every value returned by a syscall is passed to this function. It returns 0 + * if the return value is ok or -1 if there was an error. + * If the macro `VERBOSE` is defined, an appropriate message is printed to + * STDERR. + * + * @param return_value A return value from a syscall. + * + * @retval 0 The syscall was successful. + * @retval -1 There was an error. + */ +static inline signed int check_error(int return_value) { +#ifdef VERBOSE + const char *errbuf; +#endif + if (return_value < 0) { +#ifdef VERBOSE + errbuf = strerror(errno); + debug_write(errbuf); +#endif + return -1; + } + + return 0; +} + +/** + * @brief Create and connect a new TCP/IP socket + * + * This function returns a working client TCP/IP socket. + * + * @param host The host the socket will be connected to (everything resolvable, + * e.g. "::1", "8.8.8.8", "example.com") + * @param service The host's port, either numeric or as service name ("http"). + * @param proto_osi3 `LIBSOCKET_IPv4` or `LIBSOCKET_IPv6`. + * @param flags Flags to be passed to `socket(2)`. Most flags are Linux-only! + * + * @return A valid socket file descriptor. + */ +int create_inet_stream_socket(const char *host, const char *service, + char proto_osi3, int flags) { + int sfd, return_value; + struct addrinfo hint, *result, *result_check; +#ifdef VERBOSE + const char *errstring; +#endif + + if (host == NULL || service == NULL) return -1; + + memset(&hint, 0, sizeof hint); + + // set address family + switch (proto_osi3) { + case LIBSOCKET_IPv4: + hint.ai_family = AF_INET; + break; + case LIBSOCKET_IPv6: + hint.ai_family = AF_INET6; + break; + case LIBSOCKET_BOTH: + hint.ai_family = AF_UNSPEC; + break; + default: + return -1; + } + + // Transport protocol is TCP + hint.ai_socktype = SOCK_STREAM; + + if (0 != (return_value = getaddrinfo(host, service, &hint, &result))) { +#ifdef VERBOSE + errstring = gai_strerror(return_value); + debug_write(errstring); +#endif + return -1; + } + + // As described in "The Linux Programming Interface", Michael Kerrisk 2010, + // chapter 59.11 (p. 1220ff) + + for (result_check = result; result_check != NULL; + result_check = result_check->ai_next) // go through the linked list of + // struct addrinfo elements + { + sfd = socket(result_check->ai_family, result_check->ai_socktype | flags, + result_check->ai_protocol); + + if (sfd < 0) // Error!!! + continue; + + int CON_RES = connect(sfd, result_check->ai_addr, + result_check->ai_addrlen); + if ((CON_RES != -1) || (CON_RES == -1 && (flags |= SOCK_NONBLOCK) && ((errno == EINPROGRESS) || (errno == EALREADY) || (errno == EINTR)))) // connected without error, or, connected with errno being one of these important states + break; + + close(sfd); + } + + // We do now have a working socket STREAM connection to our target + + if (result_check == NULL) // Have we? + { +#ifdef VERBOSE + debug_write( + "create_inet_stream_socket: Could not connect to any address!\n"); +#endif + int errno_saved = errno; + close(sfd); + errno = errno_saved; + freeaddrinfo(result); + return -1; + } + // Yes :) + + freeaddrinfo(result); + + return sfd; +} + +/** + * @brief Creates a new UDP/IP socket + * + * Returns an integer describing a DGRAM (UDP) socket. The socket is + * automatically bound to some port. + * + * @param proto_osi3 is LIBSOCKET_IPv4 (AF_INET) or LIBSOCKET_IPv6 (AF_INET6). + * @param flags may be the flags specified in socket(2), i.e. SOCK_NONBLOCK + * and/or SOCK_CLOEXEC. More than one flags may be ORed. This argument is only + * sensible on Linux >= 2.6.27! + * + * @return The socket file descriptor number, on error -1. + * + * To send and receive data with this socket use the functions explained below, + * sendto_inet_dgram_socket() and recvfrom_inet_dgram_socket(). + */ +int create_inet_dgram_socket(char proto_osi3, int flags) { + int sfd; + + if (proto_osi3 != LIBSOCKET_IPv4 && proto_osi3 != LIBSOCKET_IPv6) { +#ifdef VERBOSE + debug_write( + "create_inet_dgram_socket: osi3 argument invalid for DGRAM " + "sockets\n"); +#endif + return -1; + } + + switch (proto_osi3) { + case LIBSOCKET_IPv4: + sfd = socket(AF_INET, SOCK_DGRAM | flags, 0); + break; + case LIBSOCKET_IPv6: + sfd = socket(AF_INET6, SOCK_DGRAM | flags, 0); + break; + default: + return -1; + } + + if (-1 == check_error(sfd)) return -1; + + return sfd; +} + +/** + * @brief This function is the equivalent to `sendto(2)` + * + * @param sfd is the *Socket File Descriptor* (every socket file descriptor + * argument in libsocket is called sfd) which you got from + * create_inet_dgram_socket(). *The usage with STREAM sockets is not recommended + * and the result is undefined!* + * @param buf is a pointer to some data. + * @param size is the length of the buffer to which buf points. + * @param host is the host to which we want to send the data. It's a string so + * you may specify everything what's resolved by getaddrinfo(), i.e. an IP + * ("193.21.34.21") or a hostname ("example.net"). + * @param service is the port on the remote host. Like in host, you may specify + * the port either as number ("123") or as service string ("ntp", "http", + * "gopher"). + * @param sendto_flags is available on all platforms. The value given here goes + * directly to the internal sendto() call. The flags which may be specified + * differ between the platforms. + * + * If it is not possible to send data at the moment, this call blocks excepted + * you specified SOCK_NONBLOCK when creating the socket. + * + * @retval n *n* bytes of data could be sent. + * @retval -1 Error. + */ +ssize_t sendto_inet_dgram_socket(int sfd, const void *buf, size_t size, + const char *host, const char *service, + int sendto_flags) { + struct sockaddr_storage oldsock; + struct addrinfo *result, *result_check, hint; + socklen_t oldsocklen = sizeof(struct sockaddr_storage); + int return_value; +#ifdef VERBOSE + const char *errstring; +#endif + + if (sfd < 0) return -1; + + if (buf == NULL) return -1; + + if (size == 0) return 0; + + if (host == NULL || service == NULL) return -1; + + if (-1 == check_error(getsockname(sfd, (struct sockaddr *)&oldsock, + (socklen_t *)&oldsocklen))) + return -1; + + memset(&hint, 0, sizeof(struct addrinfo)); + + /* + * This works for Linux > 2.6.32 + socklen_t dom_len = sizeof(hint.ai_family); + getsockopt(sfd,SOL_SOCKET,SO_DOMAIN,&hint.ai_family,&dom_len); + */ + hint.ai_family = oldsock.ss_family; + hint.ai_socktype = SOCK_DGRAM; + + if (0 != (return_value = getaddrinfo(host, service, &hint, &result))) { +#ifdef VERBOSE + errstring = gai_strerror(return_value); + debug_write(errstring); +#endif + return -1; + } + + for (result_check = result; result_check != NULL; + result_check = result_check->ai_next) // go through the linked list of + // struct addrinfo elements + { + if (-1 != (return_value = sendto( + sfd, buf, size, sendto_flags, result_check->ai_addr, + result_check->ai_addrlen))) // connected without error + { + break; // Exit loop if send operation was successful + } else { + check_error(return_value); + } + } + + freeaddrinfo(result); + + return return_value; +} + +/** + * @brief Receive data from a UDP/IP socket + * + * Receives data like `recvfrom(2)`. Pointers may be `NULL`, then the + * information (e.g. the source port) is lost (you may use NULL pointers if + * you're not interested in some information) + * + * @param sfd The socket file descriptor. + * @param buffer Where the data will be written + * @param size The size of `buffer` + * @param src_host Where the sending host's name/IP will be stored + * @param src_host_len `src_host`'s length + * @param src_service Where the port on remote side will be written to + * @param src_service_len `src_service`'s length + * @param recvfrom_flags Flags for `recvfrom(2)` + * @param numeric `LIBSOCKET_NUMERIC` if you want the names to remain + * unresolved. + * + * @retval n *n* bytes of data were received. + * @retval 0 Peer sent EOF. + * @retval <0 An error occurred. + */ +ssize_t recvfrom_inet_dgram_socket(int sfd, void *buffer, size_t size, + char *src_host, size_t src_host_len, + char *src_service, size_t src_service_len, + int recvfrom_flags, int numeric) { + struct sockaddr_storage client; + +#ifdef _TRADITIONAL_RDNS + struct sockaddr_storage oldsockaddr; + socklen_t oldsockaddrlen = sizeof(struct sockaddr_storage); + struct hostent *he; + void *addrptr; + size_t addrlen; + uint16_t sport = 0; +#endif + + ssize_t bytes; + +#ifndef _TRADITIONAL_RDNS + int retval; +#endif + +#ifdef VERBOSE + const char *errstr; +#endif + if (sfd < 0) return -1; + + if (buffer == NULL || size == 0) + return -1; + else + memset(buffer, 0, size); + + if (src_host) memset(src_host, 0, src_host_len); + if (src_service) memset(src_service, 0, src_service_len); + + socklen_t stor_addrlen = sizeof(struct sockaddr_storage); + + if (-1 == check_error(bytes = recvfrom(sfd, buffer, size, recvfrom_flags, + (struct sockaddr *)&client, + &stor_addrlen))) + return -1; + + if (src_host_len > 0 || + src_service_len > + 0) // If one of the things is wanted. If you give a null pointer + // with a positive _len parameter, you won't get the address. + { + if (numeric == LIBSOCKET_NUMERIC) { + numeric = NI_NUMERICHOST | NI_NUMERICSERV; + } + + // getnameinfo() doesn't work on FreeBSD (here) +#ifndef _TRADITIONAL_RDNS + if (0 != + (retval = getnameinfo( + (struct sockaddr *)&client, sizeof(struct sockaddr_storage), + src_host, src_host_len, src_service, src_service_len, + numeric))) // Write information to the provided memory + { +#ifdef VERBOSE + errstr = gai_strerror(retval); + debug_write(errstr); +#endif + return -1; + } +#endif + + // so use traditional methods +#ifdef _TRADITIONAL_RDNS + if (-1 == check_error(getsockname(sfd, (struct sockaddr *)&oldsockaddr, + &oldsockaddrlen))) + return -1; + + if (oldsockaddrlen > + sizeof(struct sockaddr_storage)) // If getsockname truncated the + // struct + return -1; + + if (oldsockaddr.ss_family == AF_INET) { + addrptr = &(((struct sockaddr_in *)&client)->sin_addr); + addrlen = sizeof(struct in_addr); + sport = ntohs(((struct sockaddr_in *)&client)->sin_port); + } else if (oldsockaddr.ss_family == AF_INET6) { + addrptr = &(((struct sockaddr_in6 *)&client)->sin6_addr); + addrlen = sizeof(struct in6_addr); + sport = ntohs(((struct sockaddr_in6 *)&client)->sin6_port); + } else { +#ifdef VERBOSE + debug_write("recvfrom_inet_dgram_socket: Unknown address family"); +#endif + } + + if (NULL == + (he = gethostbyaddr(addrptr, addrlen, oldsockaddr.ss_family))) { + check_error(-1); + return -1; + } + + strncpy(src_host, he->h_name, src_host_len); + snprintf(src_service, src_service_len, "%u", sport); +#endif + } + + return bytes; +} + +/** + * @brief Connect a UDP socket. + * + * If a datagram socket is connected, all data written to it (using `write(2)`) + * is sent to the peer connected to and all data `read(2)` from it is data sent + * by the peer. Usually used by clients only. + * + * @param sfd The socket file descriptor + * @param host The host to connect to + * @param service The port/service specifier + * + * @retval 0 Success + * @retval -1 Error. + */ +int connect_inet_dgram_socket(int sfd, const char *host, const char *service) { + struct addrinfo *result, *result_check, hint; + struct sockaddr_storage oldsockaddr; + socklen_t oldsockaddrlen = sizeof(struct sockaddr_storage); + int return_value; +#ifdef VERBOSE + const char *errstring; +#endif + + if (sfd < 0) return -1; + + if (host == NULL) { + // This does not work on FreeBSD systems. We pretend to disconnect the + // socket although we don't do so. This is not very severe for the + // application +#ifndef __FreeBSD__ + struct sockaddr deconnect; + memset(&deconnect, 0, sizeof(struct sockaddr)); + + deconnect.sa_family = AF_UNSPEC; + + if (check_error(connect(sfd, &deconnect, sizeof(struct sockaddr)))) + return -1; +#endif + return 0; + } + + if (-1 == check_error(getsockname(sfd, (struct sockaddr *)&oldsockaddr, + &oldsockaddrlen))) + return -1; + + if (oldsockaddrlen > + sizeof(struct sockaddr_storage)) // If getsockname truncated the struct + return -1; + + memset(&hint, 0, sizeof(struct addrinfo)); + + hint.ai_family = ((struct sockaddr_in *)&oldsockaddr) + ->sin_family; // AF_INET or AF_INET6 - offset is same + // at sockaddr_in and sockaddr_in6 + hint.ai_socktype = SOCK_DGRAM; + + if (0 != (return_value = getaddrinfo(host, service, &hint, &result))) { +#ifdef VERBOSE + errstring = gai_strerror(return_value); + debug_write(errstring); +#endif + return -1; + } + + // As described in "The Linux Programming Interface", Michael Kerrisk 2010, + // chapter 59.11 (p. 1220ff) + + for (result_check = result; result_check != NULL; + result_check = result_check->ai_next) // go through the linked list of + // struct addrinfo elements + { + if (-1 != (return_value = connect( + sfd, result_check->ai_addr, + result_check->ai_addrlen))) // connected without error + { + break; + } else { + check_error(return_value); + } + } + + // We do now have a working (updated) socket connection to our target + + if (result_check == NULL) // or not? + { +#ifdef VERBOSE + debug_write( + "connect_inet_dgram_socket: Could not connect to any address!\n"); +#endif + freeaddrinfo(result); + return -1; + } + + freeaddrinfo(result); + + return 0; +} + +/** + * @brief Close a socket. + * + * This function closes a socket. You may also use `close(2)`. + * + * @param sfd The file descriptor + * + * @retval 0 Closed socket successfully + * @retval -1 Socket was already closed (other errors are very unlikely to + * occur) + */ +int destroy_inet_socket(int sfd) { + if (sfd < 0) return -1; + + if (-1 == check_error(close(sfd))) return -1; + + return 0; +} + +/** + * @brief Perform a `shutdown(2)` call on a socket + * + * If you're done with writing or reading from a socket + * you may signalize this to the OS and/or the peer. For + * example, shutting down a socket for writing sends + * the peer an EOF signal. + * + * @param sfd The socket + * @param method `LIBSOCKET_READ` or `LIBSOCKET_WRITE` or the combination via + * `|` + * + * @retval 0 Everything's fine. + * @retval -1 Something went wrong, e.g. the socket was closed, the file + * descriptor is invalid etc. + */ +int shutdown_inet_stream_socket(int sfd, int method) { + if (sfd < 0) return -1; + + if ((method != LIBSOCKET_READ) && (method != LIBSOCKET_WRITE) && + (method != (LIBSOCKET_READ | LIBSOCKET_WRITE))) + return -1; + + if (method & LIBSOCKET_READ) // READ is set (0001 && 0001 => 0001 => true) + { + if (-1 == check_error(shutdown(sfd, SHUT_RD))) return -1; + } + + if (method & + LIBSOCKET_WRITE) // WRITE is set (0010 && 0010 => 0010 => true) + { + if (-1 == check_error(shutdown(sfd, SHUT_WR))) return -1; + } + + return 0; +} + +/* + * Server part + * + */ + +/** + * @brief Create a TCP or UDP server socket + * + * To accept connections from clients via TCP or receive datagrams via UDP, you + * need to create a server socket. This function creates such a socket and + * `bind(2)`s it to the specified address. If `proto_osi4` is `LIBSOCKET_TCP`, + * `listen(2)` is called, too. + * + * @param bind_addr Address to bind to. If you want to bind to every address use + * "0.0.0.0" or "::" (IPv6 wildcard) + * @param bind_port The port to bind to. If you write a webserver, this will be + * "http" or "80" or "https" or "443". + * @param proto_osi4 Either `LIBSOCKET_TCP` or `LIBSOCKET_UDP`. Server sockets + * in TCP and UDP differ only in that TCP sockets need a call to `listen(2)` + * @param proto_osi3 Either `LIBSOCKET_IPv4`, `LIBSOCKET_IPv6` or + * `LIBSOCKET_BOTH`; latter means that the DNS resolver should decide. + * @param flags The `flags` argument is passed ORed to the `type` argument of + * `socket(2)`; everything other than 0 does not make sense on other OSes than + * Linux. + * + * @retval >0 A working passive socket. Call `accept_inet_stream_socket()` next. + * @retval <0 Something went wrong; for example, the addresses where garbage or + * the port was not free. + */ +// Bind address Port TCP/UDP IPv4/6 +int create_inet_server_socket(const char *bind_addr, const char *bind_port, + char proto_osi4, char proto_osi3, int flags) { + int sfd, domain, type, retval; + struct addrinfo *result, *result_check, hints; +#ifdef VERBOSE + const char *errstr; +#endif + + // if ( flags != SOCK_NONBLOCK && flags != SOCK_CLOEXEC && flags != + // (SOCK_CLOEXEC|SOCK_NONBLOCK) && flags != 0 ) return -1; + + if (bind_addr == NULL || bind_port == NULL) return -1; + + switch (proto_osi4) { + case LIBSOCKET_TCP: + type = SOCK_STREAM; + break; + case LIBSOCKET_UDP: + type = SOCK_DGRAM; + break; + default: + return -1; + } + switch (proto_osi3) { + case LIBSOCKET_IPv4: + domain = AF_INET; + break; + case LIBSOCKET_IPv6: + domain = AF_INET6; + break; + case LIBSOCKET_BOTH: + domain = AF_UNSPEC; + break; + default: + return -1; + } + + memset(&hints, 0, sizeof(struct addrinfo)); + + hints.ai_socktype = type; + hints.ai_family = domain; + hints.ai_flags = AI_PASSIVE; + + if (0 != (retval = getaddrinfo(bind_addr, bind_port, &hints, &result))) { +#ifdef VERBOSE + errstr = gai_strerror(retval); + debug_write(errstr); +#endif + return -1; + } + + // As described in "The Linux Programming Interface", Michael Kerrisk 2010, + // chapter 59.11 (p. 1220ff) + for (result_check = result; result_check != NULL; + result_check = result_check->ai_next) // go through the linked list of + // struct addrinfo elements + { + sfd = socket(result_check->ai_family, result_check->ai_socktype | flags, + result_check->ai_protocol); + + if (sfd < 0) // Error at socket()!!! + continue; + + retval = bind(sfd, result_check->ai_addr, + (socklen_t)result_check->ai_addrlen); + + if (retval != 0) // Error at bind()!!! + { + close(sfd); + continue; + } + + if (type == LIBSOCKET_TCP) retval = listen(sfd, LIBSOCKET_BACKLOG); + + if (retval == 0) // If we came until here, there wasn't an error + // anywhere. It is safe to cancel the loop here + break; + else + close(sfd); + } + + if (result_check == NULL) { +#ifdef VERBOSE + debug_write( + "create_inet_server_socket: Could not bind to any address!\n"); +#endif + freeaddrinfo(result); + return -1; + } + + // We do now have a working socket on which we may call accept() + + freeaddrinfo(result); + + return sfd; +} + +/** + * @brief Accept a connection attempt on a server socket. + * + * This function accepts an incoming connection on a server socket. + * + * (the `src_*` arguments may be `NULL`, in which case the address is not + * stored) + * + * @param sfd The server socket + * @param src_host Buffer where the client's address is copied to + * @param src_host_len `src_host`'s length. If the hostname is longer than this, + * it is truncated. + * @param src_service Buffer in which the client's port is stored + * @param src_service_len Its size. If shorter than the hostname it gets + * truncated. + * @param flags May be `LIBSOCKET_NUMERIC`; then there is no rDNS lookup and the + * IP and port number are stored as-is. + * @param accept_flags Flags for `accept4(2)` (which is only used on Linux) + * + * @retval >0 A socket file descriptor which can be used to talk to the client + * @retval <0 Error. + */ +// Socket Src string Src str len Src service +// Src service len NUMERIC? +int accept_inet_stream_socket(int sfd, char *src_host, size_t src_host_len, + char *src_service, size_t src_service_len, + int flags, int accept_flags) { + struct sockaddr_storage client_info; + int client_sfd; + +#ifndef _TRADITIONAL_RDNS + int retval; +#endif + +#ifdef _TRADITIONAL_RDNS + struct sockaddr_storage oldsockaddr; + socklen_t oldsockaddrlen = sizeof(struct sockaddr_storage); + struct hostent *he; + void *addrptr; + size_t in_addrlen; + uint16_t sport = 0; +#endif + +#ifdef VERBOSE + const char *errstr; +#endif + socklen_t addrlen = sizeof(struct sockaddr_storage); + + // Portable behavior +#if LIBSOCKET_LINUX + if (-1 == + check_error((client_sfd = accept4(sfd, (struct sockaddr *)&client_info, + &addrlen, accept_flags)))) // blocks + return -1; +#else + if (-1 == + check_error((client_sfd = accept(sfd, (struct sockaddr *)&client_info, + &addrlen)))) // blocks + return -1; +#endif + + if (src_host_len > 0 || + src_service_len > + 0) // If one of the things is wanted. If you give a null pointer + // with a positive _len parameter, you won't get the address. + { + if (flags == LIBSOCKET_NUMERIC) { + flags = NI_NUMERICHOST | NI_NUMERICSERV; + } else { + flags = 0; // To prevent errors: Unknown flags are ignored + } + +#ifndef _TRADITIONAL_RDNS + if (0 != (retval = getnameinfo( + (struct sockaddr *)&client_info, + sizeof(struct sockaddr_storage), src_host, src_host_len, + src_service, src_service_len, + flags))) // Write information to the provided memory + { +#ifdef VERBOSE + errstr = gai_strerror(retval); + debug_write(errstr); +#endif + } +#endif + +#ifdef _TRADITIONAL_RDNS + if (-1 == check_error(getsockname(sfd, (struct sockaddr *)&oldsockaddr, + &oldsockaddrlen))) + return -1; + + if (oldsockaddrlen > + sizeof(struct sockaddr_storage)) // If getsockname truncated the + // struct + return -1; + + if (oldsockaddr.ss_family == AF_INET) { + addrptr = &(((struct sockaddr_in *)&client_info)->sin_addr); + in_addrlen = sizeof(struct in_addr); + sport = ntohs(((struct sockaddr_in *)&client_info)->sin_port); + } else if (oldsockaddr.ss_family == AF_INET6) { + addrptr = &(((struct sockaddr_in6 *)&client_info)->sin6_addr); + in_addrlen = sizeof(struct in6_addr); + sport = ntohs(((struct sockaddr_in6 *)&client_info)->sin6_port); + } else { +#ifdef VERBOSE + debug_write("accept_inet_stream_socket: Unknown address family"); +#endif + return -1; + } + + if (NULL == + (he = gethostbyaddr(addrptr, in_addrlen, oldsockaddr.ss_family))) { + check_error(-1); + // Don't return with error on name resolution failure + } + + strncpy(src_host, he->h_name, src_host_len); + snprintf(src_service, src_service_len, "%u", sport); +#endif + } + + return client_sfd; +} + +/** + * @brief Look up which address families a host supports. + * + * If you want to send a datagram to a host but you don't know + * if it supports IPv4 or IPv6, use this function. It returns + * the address family returned by a DNS lookup. On most systems IPv6 + * is the preferred address family. + * + * @param hostname The hostname of the host you want to look up. + * + * @retval LIBSOCKET_IPv4 Host supports only IPv4 + * @retval LIBSOCKET_IPv6 Host supports IPv6 (usually it supports IPv4 then, + * too) + * @retval <0 Error. + */ +int get_address_family(const char *hostname) { + int return_value; + struct addrinfo hint, *result; +#ifdef VERBOSE + const char *errstring; +#endif + int af; + + if (hostname == NULL) return -1; + + memset(&hint, 0, sizeof hint); + + hint.ai_family = AF_UNSPEC; + + if (0 != (return_value = getaddrinfo(hostname, "0", &hint, &result))) { +#ifdef VERBOSE + errstring = gai_strerror(return_value); + debug_write(errstring); +#endif + return -1; + } + + if (result == NULL) return -1; + + if (result->ai_family == AF_INET) { + af = LIBSOCKET_IPv4; + } else if (result->ai_family == AF_INET6) { + af = LIBSOCKET_IPv6; + } else { + af = -1; + } + + return af; +} + +/** + * @brief Create a datagram socket and join to the multicast group `address`. + * + * This function creates a multicast socket bound to `address`. The only option + * set is the `IP_MULTICAST_IF` (`IPV6_MULTICAST_IF`) option to avoid an + * explicit routing entry for the multicast address. + * + * An option you want to set very likely is `IP_MULTICAST_LOOP`. Refer to + * `ip(7)` respectively `ipv6(7)` for `setsockopt()` options; for example: + * + * int c = 0; + * setsockopt(sfd,IPPROTO_IP,IP_MULTICAST_LOOP,&c,4); + * + * The group address and port is also used as arguments to `bind(2)`. After + * creating this socket, you may use the usual I/O functions on it, i.e. + * sendto_inet_dgram_socket and recvfrom_inet_dgram_socket. + * + * @param group Group address. This address is also used to bind the socket + * @param port Multicast port. + * @param local For IPv4 multicast groups: The address of the interface to be + * used. Ignored for IPv6, NULL for kernel's choice + * + * @retval <0 Error (Check errno or use `LIBSOCKET_VERBOSE`) + * @retval >=0 A valid file descriptor. + * + */ +#ifdef LIBSOCKET_LINUX +int create_multicast_socket(const char *group, const char *port, + const char *if_name) { + int sfd, return_value; + struct sockaddr maddr, localif; + struct addrinfo hints, *result; + struct ip_mreqn mreq4; + struct ipv6_mreq mreq6; + struct in_addr any; + struct ifreq interface; + + memset(&maddr, 0, sizeof(maddr)); + memset(&localif, 0, sizeof(localif)); + memset(&mreq4, 0, sizeof(mreq4)); + memset(&mreq6, 0, sizeof(mreq6)); + memset(&hints, 0, sizeof(hints)); + memset(&interface, 0, sizeof(interface)); + + if (-1 == check_error(sfd = create_inet_server_socket( + group, port, LIBSOCKET_UDP, LIBSOCKET_BOTH, 0))) { + return -1; + } + + hints.ai_socktype = SOCK_DGRAM; + hints.ai_family = AF_UNSPEC; + + if (0 != (return_value = getaddrinfo(group, port, &hints, &result))) { + int errno_saved = errno; +#ifdef VERBOSE + const char *errstring = gai_strerror(return_value); + debug_write(errstring); +#endif + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + + return -1; + } + + if (result->ai_family == AF_INET) { + // Result is IPv4 address. + mreq4.imr_multiaddr = ((struct sockaddr_in *)result->ai_addr)->sin_addr; + + if (if_name == NULL) { + mreq4.imr_ifindex = 0; + any.s_addr = INADDR_ANY; + mreq4.imr_address = any; + } else { + memcpy(interface.ifr_name, if_name, + strlen(if_name) > IFNAMSIZ ? IFNAMSIZ : strlen(if_name)); + + if (-1 == check_error(ioctl(sfd, SIOCGIFINDEX, &interface))) { + int errno_saved = errno; + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + return -1; + } + + mreq4.imr_ifindex = interface.ifr_ifindex; + } + + if (-1 == check_error(setsockopt(sfd, IPPROTO_IP, IP_ADD_MEMBERSHIP, + &mreq4, sizeof(struct ip_mreqn)))) { + int errno_saved = errno; + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + return -1; + } + if (-1 == check_error(setsockopt(sfd, IPPROTO_IP, IP_MULTICAST_IF, + &mreq4, sizeof(struct ip_mreqn)))) { + int errno_saved = errno; + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + return -1; + } + + // Setup finished. + // + } else if (result->ai_family == AF_INET6) { + mreq6.ipv6mr_multiaddr = + ((struct sockaddr_in6 *)result->ai_addr)->sin6_addr; + mreq6.ipv6mr_interface = 0; + + if (if_name == NULL) + mreq6.ipv6mr_interface = 0; + else { + memcpy(interface.ifr_name, if_name, + strlen(if_name) > IFNAMSIZ ? IFNAMSIZ : strlen(if_name)); + + if (-1 == check_error(ioctl(sfd, SIOCGIFINDEX, &interface))) { + int errno_saved = errno; + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + return -1; + } + + mreq6.ipv6mr_interface = interface.ifr_ifindex; + } + + if (-1 == check_error(setsockopt(sfd, IPPROTO_IPV6, IPV6_ADD_MEMBERSHIP, + &mreq6, sizeof(struct ipv6_mreq)))) { + int errno_saved = errno; + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + return -1; + } + if (-1 == check_error(setsockopt(sfd, IPPROTO_IPV6, IPV6_MULTICAST_IF, + &mreq6.ipv6mr_interface, + sizeof(mreq6.ipv6mr_interface)))) { + int errno_saved = errno; + close(sfd); + freeaddrinfo(result); + errno = errno_saved; + return -1; + } + } + + freeaddrinfo(result); + return sfd; +} + +#endif + +/** + * @} + */ + +#undef debug_write diff --git a/src/socket/libinetsocket.h b/src/socket/libinetsocket.h new file mode 100644 index 00000000..71375a2e --- /dev/null +++ b/src/socket/libinetsocket.h @@ -0,0 +1,97 @@ +#ifndef LIBSOCKET_LIBINETSOCKET_H_C1A9FFEDF5E94B2FB010A0FAA0E92A2F +#define LIBSOCKET_LIBINETSOCKET_H_C1A9FFEDF5E94B2FB010A0FAA0E92A2F +/** + * @file libinetsocket.h + * + * @brief A static copy of the libsocket library header + * @authors T.Tian + * + * Original license information: + */ +/* + The committers of the libsocket project, all rights reserved + (c) 2012, dermesser + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above + copyright notice, this list of conditions and the following disclaimer. + 2. Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS “AS IS” AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY + DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +/* Headers (e.g. for flags) */ +#include +#include + +/* Macro definitions */ + +#define LIBSOCKET_TCP 1 +#define LIBSOCKET_UDP 2 + +#define LIBSOCKET_IPv4 3 +#define LIBSOCKET_IPv6 4 + +#define LIBSOCKET_BOTH \ + 5 /* let the resolver/library choose (TCP/UDP or IPv4/6) */ + +#define LIBSOCKET_READ 1 +#define LIBSOCKET_WRITE 2 + +#define LIBSOCKET_NUMERIC 1 + +#ifdef __cplusplus +#ifdef MIXED +extern "C" { +#endif +#endif + +extern int create_inet_stream_socket(const char* host, const char* service, + char proto_osi3, int flags); +extern int create_inet_dgram_socket(char proto_osi3, int flags); +extern ssize_t sendto_inet_dgram_socket(int sfd, const void* buf, size_t size, + const char* host, const char* service, + int sendto_flags); +extern ssize_t recvfrom_inet_dgram_socket(int sfd, void* buffer, size_t size, + char* src_host, size_t src_host_len, + char* src_service, + size_t src_service_len, + int recvfrom_flags, int numeric); +extern int connect_inet_dgram_socket(int sfd, const char* host, + const char* service); +extern int destroy_inet_socket(int sfd); +extern int shutdown_inet_stream_socket(int sfd, int method); +extern int create_inet_server_socket(const char* bind_addr, + const char* bind_port, char proto_osi4, + char proto_osi3, int flags); +extern int accept_inet_stream_socket(int sfd, char* src_host, + size_t src_host_len, char* src_service, + size_t src_service_len, int flags, + int accept_flags); +extern int get_address_family(const char* hostname); + +#ifdef __linux__ +extern int create_multicast_socket(const char* group, const char* port, + const char* local); +#endif + +#ifdef __cplusplus +#ifdef MIXED +} +#endif +#endif + +#endif diff --git a/src/socket/libunixsocket.c b/src/socket/libunixsocket.c new file mode 100644 index 00000000..1918a83a --- /dev/null +++ b/src/socket/libunixsocket.c @@ -0,0 +1,515 @@ +#ifndef _GNU_SOURCE +#define _GNU_SOURCE // accept4() +#endif + + +#include "conf.h" + +#include +#include +#include +#include +#include +#include +#include // UNIX domain sockets +#include // read()/write() + +/** + * @addtogroup libunixsocket + * @{ + * @file libunixsocket.c + * @brief A static copy of the libsocket library's unitsocket code. + * @authors T.Tian + * + * Original license information: + * Contains all C functions to handle UNIX domain sockets. + + The committers of the libsocket project, all rights reserved + (c) 2012, dermesser + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS “AS IS” AND ANY + EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY + DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +/* + * Structure of the functions defined here: + * + * + * + * + * + */ + +//# define VERBOSE + +// Macro definitions + +#define LIBSOCKET_BACKLOG 128 + +#define LIBSOCKET_STREAM 1 +#define LIBSOCKET_DGRAM 2 + +#define LIBSOCKET_READ 1 +#define LIBSOCKET_WRITE 2 + +#define debug_write(str) \ + { \ + int verbose_errno_save = errno; \ + write(2, str, strlen(str)); \ + errno = verbose_errno_save; \ + } + +/** + * @brief Checks return value for error. + * + * Every value returned by a syscall is passed to this function. It returns 0 + * if the return value is ok or -1 if there was an error. + * If the macro `VERBOSE` is defined, an appropriate message is printed to + * STDERR. + * + * @param return_value A return value from a syscall. + * + * @retval 0 The syscall was successful. + * @retval -1 There was an error. + */ +static inline signed int check_error(int return_value) { +#ifdef VERBOSE + const char* errbuf; +#endif + if (return_value < 0) { +#ifdef VERBOSE + errbuf = strerror(errno); + debug_write(errbuf); +#endif + return -1; + } + + return 0; +} + +static int set_unix_socket_path(struct sockaddr_un* saddr, const char* path_or_name) { + memset(&saddr->sun_path, 0, sizeof(saddr->sun_path)); + if (path_or_name[0] != 0) { + strncpy(saddr->sun_path, path_or_name, sizeof(saddr->sun_path) - 1); + } else { + // Abstract socket address. + // Make sure this is not a nulled memory region. + if (path_or_name[1] == 0) { + errno = EINVAL; +#ifdef VERBOSE + debug_write("set_unix_socket_path: Socket address is too many 0s\n"); +#endif + return -1; + } + size_t max_len = 1 + strlen(path_or_name+1); + if (max_len > sizeof(saddr->sun_path) - 1) { +#ifdef VERBOSE + debug_write("set_unix_socket_path: Abstract socket address is too long\n"); +#endif + errno = ENAMETOOLONG; + return -1; + } + memcpy(saddr->sun_path, path_or_name, max_len); + } + return 0; +} + +/** + * @brief Create and connect a new UNIX STREAM socket. + * + * Creates and connects a new STREAM socket with the socket given in `path`. + * + * If the first byte in path is a null byte, this function assumes that it is an + * abstract socket address (see `unix(7)`; this is a Linux-specific feature). + * Following that byte, there must be a normal 0-terminated string. + * + * @retval >0 Success; return value is a socket file descriptor + * @retval <0 Error. + */ +int create_unix_stream_socket(const char* path, int flags) { + struct sockaddr_un saddr; + int sfd; + + if (path == NULL) return -1; + + if (-1 == check_error(sfd = socket(AF_UNIX, SOCK_STREAM | flags, 0))) + return -1; + + memset(&saddr, 0, sizeof(struct sockaddr_un)); + + if (strlen(path) > (sizeof(saddr.sun_path) - 1)) { +#ifdef VERBOSE + debug_write( + "create_unix_stream_socket: UNIX socket destination path too " + "long\n"); +#endif + return -1; + } + + saddr.sun_family = AF_UNIX; + if (-1 == check_error(set_unix_socket_path(&saddr, path))) + return -1; + size_t pathlen = strlen(saddr.sun_path); + if (pathlen == 0) // likely abstract socket address + pathlen = sizeof(saddr.sun_path); + + + if (-1 == check_error( + connect(sfd, (struct sockaddr*)&saddr, + sizeof(saddr.sun_family) + pathlen))) { + close(sfd); + return -1; + } + + return sfd; +} + +/** + * @brief Create a UNIX DGRAM socket + * + * @param bind_path If not `NULL`, bind to `bind_path`. + * @param flags Flags to pass to `socket(2)` (varies from OS to OS; look in the + * man pages) + * + * @retval >0 Success. Value is socket. + * @retval <0 Error. + */ +int create_unix_dgram_socket(const char* bind_path, int flags) { + int sfd, retval; + struct sockaddr_un saddr; + + if (-1 == check_error(sfd = socket(AF_UNIX, SOCK_DGRAM | flags, 0))) + return -1; + + memset(&saddr, 0, sizeof(struct sockaddr_un)); + + if (bind_path != NULL) { + if ((retval = unlink(bind_path)) == -1 && + errno != ENOENT) // If there's another error than "doesn't exist" + { + check_error(retval); + return -1; + } + + if (strlen(bind_path) > (sizeof(saddr.sun_path) - 1)) { +#ifdef VERBOSE + debug_write( + "create_unix_dgram_socket: UNIX socket path too long\n"); +#endif + return -1; + } + + saddr.sun_family = AF_UNIX; + if (-1 == check_error(set_unix_socket_path(&saddr, bind_path))) + return -1; + + size_t pathlen = strlen(saddr.sun_path); + if (pathlen == 0) // likely abstract socket address + pathlen = sizeof(saddr.sun_path); + + bind(sfd, (struct sockaddr*)&saddr, + sizeof(saddr.sun_family) + pathlen); + } + + return sfd; +} + +/** + * @brief Connect a datagram socket + * + * Connects a datagram socket to the specified socket so the + * stream i/o operations may be used (`read(2)/write(2)`) + * + * @param sfd The socket + * @param path The path to connect to + * + * @retval 0 Fine + * @retval <0 Not fine + */ +int connect_unix_dgram_socket(int sfd, const char* path) { + struct sockaddr_un new_addr; + struct sockaddr deconnect; + + if (sfd < 0) return -1; + + if (path == NULL) { + memset(&deconnect, 0, sizeof(struct sockaddr)); + + deconnect.sa_family = AF_UNSPEC; + + if (check_error(connect(sfd, &deconnect, sizeof(struct sockaddr)))) + return -1; + + return 0; + } + + memset(&new_addr, 0, sizeof(struct sockaddr_un)); + + new_addr.sun_family = AF_UNIX; + + if (strlen(path) > sizeof(new_addr.sun_path) - 1) { +#ifdef VERBOSE + debug_write("connect_unix_dgram_socket: Path too long\n"); +#endif + return -1; + } + + if (-1 == check_error(set_unix_socket_path(&new_addr, path))) { + return -1; + } + + size_t pathlen = strlen(new_addr.sun_path); + if (pathlen == 0) // likely abstract socket address + pathlen = sizeof(new_addr.sun_path); + + if (-1 == check_error(connect( + sfd, (struct sockaddr*)&new_addr, + sizeof(new_addr.sun_family) + pathlen))) + return -1; + + return 0; +} + +/** + * @brief Close a socket + * + * Actually, it's the same as `close(2)`. + * + * @param sfd The socket to close. + * + * @retval 0 Socket could be closed + * @retval <0 Socket was already closed. + */ +int destroy_unix_socket(int sfd) { + if (sfd < 0) return -1; + + if (-1 == check_error(close(sfd))) return -1; + + return 0; +} + +/** + * @brief Shut a socket down + * + * Shut a socket down for reading or writing. If shut down for + * reading, the program can't read anymore. If shut down for writing + * no data can be sent anymore and the peer receives EOF. + * + * @param sfd The socket + * @param method `LIBSOCKET_READ`, `LIBSOCKET_WRITE` or + * `LIBSOCKET_READ|LIBSOCKET_WRITE` + * + * @retval 0 Success + * @retval <0 Error + */ +int shutdown_unix_stream_socket(int sfd, int method) { + if (sfd < 0) return -1; + + if ((method != LIBSOCKET_READ) && (method != LIBSOCKET_WRITE) && + (method != (LIBSOCKET_READ | LIBSOCKET_WRITE))) + return -1; + + if (method & LIBSOCKET_READ) // READ is set (0001 && 0001 => 0001) + { + if (-1 == check_error(shutdown(sfd, SHUT_RD))) return -1; + } + + if (method & LIBSOCKET_WRITE) // WRITE is set (0010 && 0010 => 0010) + { + if (-1 == check_error(shutdown(sfd, SHUT_WR))) return -1; + } + + return 0; +} + +/** + * @brief Create a passive UNIX socket + * + * Creating a DGRAM server socket is the same as creating one + * using `create_unix_dgram_socket()` but with latter you may + * also not bind to anywhere. + * + * @param path Path to bind the socket to + * @param socktype `LIBSOCKET_STREAM` or `LIBSOCKET_DGRAM` + * @param flags Flags for `socket(2)`. + * + * @retval >0 Success; returned value is a file descriptor for the socket + * @retval <0 An error occurred. + */ +int create_unix_server_socket(const char* path, int socktype, int flags) { + struct sockaddr_un saddr; + int sfd, type, retval; + + if (path == NULL) return -1; + + if (strlen(path) > (sizeof(saddr.sun_path) - 1)) { +#ifdef VERBOSE + debug_write("create_unix_server_socket: Path too long\n"); +#endif + return -1; + } + + switch (socktype) { + case LIBSOCKET_STREAM: + type = SOCK_STREAM; + break; + case LIBSOCKET_DGRAM: + type = SOCK_DGRAM; + break; + default: + return -1; + } + + if (-1 == check_error(sfd = socket(AF_UNIX, type | flags, 0))) return -1; + + if ((retval = unlink(path)) == -1 && + errno != ENOENT) // If there's another error than "doesn't exist" + { + check_error(retval); + return -1; + } + + memset(&saddr, 0, sizeof(struct sockaddr_un)); + + saddr.sun_family = AF_UNIX; + + if (-1 == check_error(set_unix_socket_path(&saddr, path))) + return -1; + + size_t pathlen = strlen(saddr.sun_path); + if (pathlen == 0) // likely abstract socket address + pathlen = sizeof(saddr.sun_path); + + if (-1 == + check_error(bind(sfd, (struct sockaddr*)&saddr, + sizeof(saddr.sun_family) + pathlen))) + return -1; + + if (type == SOCK_STREAM) { + if (-1 == check_error(listen(sfd, LIBSOCKET_BACKLOG))) return -1; + } + + return sfd; +} + +/** + * @brief Accept connections on a passive UNIX socket + * + * @param sfd The server socket + * @param flags Flags for `accept4(3)` (therefore useless on any other system + * than Linux) + * + * @retval >0 Return value is a socket connected to the client + * @retval <0 Error at `accept[4]()` + */ +int accept_unix_stream_socket(int sfd, int flags) { + int cfd; + + if (sfd < 0) return -1; +#if LIBSOCKET_LINUX + if (-1 == check_error(cfd = accept4(sfd, 0, 0, flags))) return -1; +#else + if (-1 == check_error(cfd = accept(sfd, 0, 0))) return -1; +#endif + return cfd; +} + +/** + * @brief Receive datagram from another UNIX socket + * + * @param sfd The socket descriptor + * @param buf The buffer to which the data is written + * @param size its size + * @param from Place where the path of the sending socket is placed to + * @param from_size its size + * @param recvfrom_flags Flags passed to `recvfrom(2)` + * + * @retval n *n* bytes were received + * @retval <0 Error at `recvfrom(2)` + */ +ssize_t recvfrom_unix_dgram_socket(int sfd, void* buf, size_t size, char* from, + size_t from_size, int recvfrom_flags) { + int bytes; + socklen_t socksize = sizeof(struct sockaddr_un); + struct sockaddr_un saddr; + + memset(buf, 0, size); + memset(from, 0, from_size); + + if (-1 == + check_error(bytes = recvfrom(sfd, buf, size, recvfrom_flags, + (struct sockaddr*)&saddr, &socksize))) + return -1; + + if (from != NULL && from_size > 0) { + memcpy(from, saddr.sun_path, + from_size < sizeof(saddr.sun_path) ? from_size + : sizeof(saddr.sun_path)); + } + + return bytes; +} + +/** + * @brief Send datagram to socket + * + * @param sfd Socket + * @param buf Data to be sent + * @param size The length of the buffer `buf` + * @param path Destination socket + * @param sendto_flags Flags passed to `sendto(2)` + * + * @retval n *n* bytes were sent + * @retval <0 Error at `sendto(2)`. + */ +ssize_t sendto_unix_dgram_socket(int sfd, const void* buf, size_t size, + const char* path, int sendto_flags) { + int bytes; + struct sockaddr_un saddr; + + if (strlen(path) > sizeof(saddr.sun_path) - 1) { +#ifdef VERBOSE + debug_write( + "sendto_unix_dgram_socket: UNIX destination socket path too " + "long\n"); +#endif + return -1; + } + + memset(&saddr, 0, sizeof(struct sockaddr_un)); + + saddr.sun_family = AF_UNIX; + if (-1 == check_error(set_unix_socket_path(&saddr, path))) + return -1; + + if (-1 == check_error(bytes = sendto(sfd, buf, size, sendto_flags, + (struct sockaddr*)&saddr, + sizeof(struct sockaddr_un)))) + return -1; + + return bytes; +} + +/** + * @} + */ + +#undef debug_write diff --git a/src/socket/libunixsocket.h b/src/socket/libunixsocket.h new file mode 100644 index 00000000..22df48f5 --- /dev/null +++ b/src/socket/libunixsocket.h @@ -0,0 +1,73 @@ +#ifndef LIBSOCKET_LIBUNIXSOCKET_H_61CF2FC7034E4AD982DA08144D578572 +#define LIBSOCKET_LIBUNIXSOCKET_H_61CF2FC7034E4AD982DA08144D578572 +/** + * @file libunixsocket.h + * @brief A static copy of the libsocket library header + * @authors T.Tian + * + * Original license information: + * Contains all libunixsocket C functions. + * + The committers of the libsocket project, all rights reserved + (c) 2012, dermesser + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above + copyright notice, this list of conditions and the following disclaimer. + 2. Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS “AS IS” AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY + DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +/* Headers (e.g. for flags) */ +#include +#include + +/* Macro definitions */ + +#define LIBSOCKET_STREAM 1 +#define LIBSOCKET_DGRAM 2 + +#define LIBSOCKET_READ 1 +#define LIBSOCKET_WRITE 2 + +#ifdef __cplusplus +#ifdef MIXED +extern "C" { +#endif +#endif + +extern int create_unix_stream_socket(const char* path, int flags); +extern int create_unix_dgram_socket(const char* bind_path, int flags); +extern int connect_unix_dgram_socket(int sfd, const char* path); +extern int destroy_unix_socket(int sfd); +extern int shutdown_unix_stream_socket(int sfd, int method); +extern int create_unix_server_socket(const char* path, int socktype, int flags); +extern int accept_unix_stream_socket(int sfd, int flags); +extern ssize_t recvfrom_unix_dgram_socket(int sfd, void* buf, size_t size, + char* from, size_t from_size, + int recvfrom_flags); +extern ssize_t sendto_unix_dgram_socket(int sfd, const void* buf, size_t size, + const char* path, int sendto_flags); + +#ifdef __cplusplus +#ifdef MIXED +} +#endif +#endif + +#endif From 8cfeeeb06f7355b6579913f0b429a16e1557e510 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 17:11:31 +0800 Subject: [PATCH 02/14] revert diff checkers --- src/md.c | 265 +------------------------------------------------------ 1 file changed, 1 insertion(+), 264 deletions(-) diff --git a/src/md.c b/src/md.c index 7f9bd2ef..028c0281 100644 --- a/src/md.c +++ b/src/md.c @@ -712,7 +712,6 @@ void NVK_G(SPARC_OBJ *pSPARC) { // Charge extrapolation (for better rho_guess) elecDensExtrapolation(pSPARC); - // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem @@ -854,20 +853,14 @@ void NPT_NH (SPARC_OBJ *pSPARC) { MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&pSPARC->pres_i, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); -<<<<<<< variant A ->>>>>>> variant B -======= end if ((pSPARC->MDCount != 1) || (pSPARC->RestartFlag == 1)) { // Update velocity of particles in the second half timestep AccelVelocityParticle(pSPARC); // Update velocity of virtual thermo and baro variables in the second half timestep IsoPress(pSPARC); } -<<<<<<< variant A ->>>>>>> variant B -======= end // Update velocity of virtual thermo and baro variables in the first half timestep IsoPress(pSPARC); // Update velocity of particles in the first half timestep @@ -876,26 +869,14 @@ void NPT_NH (SPARC_OBJ *pSPARC) { PositionParticleCell(pSPARC); // Reinitialize mesh size and related variables after changing size of cell reinitialize_mesh_NPT(pSPARC); -<<<<<<< variant A ->>>>>>> variant B -======= end // Charge extrapolation (for better rho_guess) elecDensExtrapolation(pSPARC); // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); -<<<<<<< variant A - // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF - if (pSPARC->mlff_flag == 0){ - Calculate_electronicGroundState(pSPARC); - } else { - MLFF_call(pSPARC); - } ->>>>>>> variant B // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem Calculate_Properties(pSPARC); //Calculate_electronicGroundState(pSPARC); -======= end pSPARC->elecgs_Count++; #ifdef DEBUG // Calculate Hamiltonian of the system. @@ -904,10 +885,7 @@ void NPT_NH (SPARC_OBJ *pSPARC) { printf("\nend NPT timestep %d\n", pSPARC->MDCount + 1); } #endif -<<<<<<< variant A ->>>>>>> variant B -======= end } /* @@ -1139,11 +1117,7 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) double t1, t2; #endif -<<<<<<< variant A int p, i; ->>>>>>> variant B - int p; -======= end // scaling factor double scal = pSPARC->scale; // the ratio between length #ifdef DEBUG @@ -1156,13 +1130,9 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->delta_y = pSPARC->range_y/(pSPARC->numIntervals_y); pSPARC->delta_z = pSPARC->range_z/(pSPARC->numIntervals_z); -<<<<<<< variant A pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z * pSPARC->Jacbdet; ->>>>>>> variant B -======= end -<<<<<<< variant A #ifdef DEBUG if (!rank) { // printf("Volume: %12.6f\n", vol); @@ -1173,18 +1143,8 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) if (i%3==2 && i>0) printf("\n"); } printf("\n"); ->>>>>>> variant B - if (pSPARC->cell_typ <= 20) { - pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z * pSPARC->Jacbdet; - } else if (pSPARC->cell_typ > 20 && pSPARC->cell_typ <= 30) { - pSPARC->dV = pSPARC->delta_x * pSPARC->delta_y * pSPARC->delta_z; - // pSPARC->twist = pSPARC->twistpercell/pSPARC->range_z; -======= end } -<<<<<<< variant A #endif ->>>>>>> variant B -======= end int FDn = pSPARC->order / 2; @@ -1199,7 +1159,6 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->D1_stencil_coeffs_z[p] = pSPARC->FDweights_D1[p] * dz_inv; } -<<<<<<< variant A // 2nd derivative weights including mesh double dx2_inv, dy2_inv, dz2_inv; dx2_inv = 1.0 / (pSPARC->delta_x * pSPARC->delta_x); @@ -1229,21 +1188,7 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->D1_stencil_coeffs_zy[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // d/dz(2*T_23 df/dy) used in d/dz(2*T_12 df/dx + 2*T_23 df/dy) } } ->>>>>>> variant B - if (pSPARC->CyclixFlag) { - LapStencil_cyclix(pSPARC); - Integration_weights_cyclix(pSPARC, pSPARC->Intgwt_kpttopo, pSPARC->DMVertices_kptcomm[0], pSPARC->Nx_d_kptcomm, pSPARC->Ny_d_kptcomm, pSPARC->Nz_d_kptcomm); - Integration_weights_cyclix(pSPARC, pSPARC->Intgwt_psi, pSPARC->DMVertices_dmcomm[0], pSPARC->Nx_d_dmcomm, pSPARC->Ny_d_dmcomm, pSPARC->Nz_d_dmcomm); - Integration_weights_cyclix(pSPARC, pSPARC->Intgwt_phi, pSPARC->DMVertices[0], pSPARC->Nx_d, pSPARC->Ny_d, pSPARC->Nz_d); - } else { - // 2nd derivative weights including mesh - double dx2_inv, dy2_inv, dz2_inv; - dx2_inv = 1.0 / (pSPARC->delta_x * pSPARC->delta_x); - dy2_inv = 1.0 / (pSPARC->delta_y * pSPARC->delta_y); - dz2_inv = 1.0 / (pSPARC->delta_z * pSPARC->delta_z); -======= end - -<<<<<<< variant A + // maximum eigenvalue of -0.5 * Lap (only with periodic boundary conditions) if(pSPARC->cell_typ == 0) { #ifdef DEBUG @@ -1260,32 +1205,7 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->MaxEigVal_mhalfLap += 2.0 * (pSPARC->D2_stencil_coeffs_x[p] * cos(M_PI*p*scal_x) + pSPARC->D2_stencil_coeffs_y[p] * cos(M_PI*p*scal_y) + pSPARC->D2_stencil_coeffs_z[p] * cos(M_PI*p*scal_z)); ->>>>>>> variant B - // Stencil coefficients for mixed derivatives - if (pSPARC->cell_typ == 0) { - for (p = 0; p < FDn + 1; p++) { - pSPARC->D2_stencil_coeffs_x[p] = pSPARC->FDweights_D2[p] * dx2_inv; - pSPARC->D2_stencil_coeffs_y[p] = pSPARC->FDweights_D2[p] * dy2_inv; - pSPARC->D2_stencil_coeffs_z[p] = pSPARC->FDweights_D2[p] * dz2_inv; - } - } else if (pSPARC->cell_typ > 10 && pSPARC->cell_typ < 20) { - for (p = 0; p < FDn + 1; p++) { - pSPARC->D2_stencil_coeffs_x[p] = pSPARC->lapcT[0] * pSPARC->FDweights_D2[p] * dx2_inv; - pSPARC->D2_stencil_coeffs_y[p] = pSPARC->lapcT[4] * pSPARC->FDweights_D2[p] * dy2_inv; - pSPARC->D2_stencil_coeffs_z[p] = pSPARC->lapcT[8] * pSPARC->FDweights_D2[p] * dz2_inv; - pSPARC->D2_stencil_coeffs_xy[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dx_inv; // 2*T_12 d/dx(df/dy) - pSPARC->D2_stencil_coeffs_xz[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dx_inv; // 2*T_13 d/dx(df/dz) - pSPARC->D2_stencil_coeffs_yz[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // 2*T_23 d/dy(df/dz) - pSPARC->D1_stencil_coeffs_xy[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dy_inv; // d/dx(2*T_12 df/dy) used in d/dx(2*T_12 df/dy + 2*T_13 df/dz) - pSPARC->D1_stencil_coeffs_yx[p] = 2 * pSPARC->lapcT[1] * pSPARC->FDweights_D1[p] * dx_inv; // d/dy(2*T_12 df/dx) used in d/dy(2*T_12 df/dx + 2*T_23 df/dz) - pSPARC->D1_stencil_coeffs_xz[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dz_inv; // d/dx(2*T_13 df/dz) used in d/dx(2*T_12 df/dy + 2*T_13 df/dz) - pSPARC->D1_stencil_coeffs_zx[p] = 2 * pSPARC->lapcT[2] * pSPARC->FDweights_D1[p] * dx_inv; // d/dz(2*T_13 df/dx) used in d/dz(2*T_13 df/dz + 2*T_23 df/dy) - pSPARC->D1_stencil_coeffs_yz[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dz_inv; // d/dy(2*T_23 df/dz) used in d/dy(2*T_12 df/dx + 2*T_23 df/dz) - pSPARC->D1_stencil_coeffs_zy[p] = 2 * pSPARC->lapcT[5] * pSPARC->FDweights_D1[p] * dy_inv; // d/dz(2*T_23 df/dy) used in d/dz(2*T_12 df/dx + 2*T_23 df/dy) - } -======= end } -<<<<<<< variant A pSPARC->MaxEigVal_mhalfLap *= -0.5; #ifdef DEBUG t2 = MPI_Wtime(); @@ -1293,10 +1213,7 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) pSPARC->MaxEigVal_mhalfLap, (t2-t1)*1e3); #endif } ->>>>>>> variant B -======= end -<<<<<<< variant A double h_eff = 0.0; if (fabs(pSPARC->delta_x - pSPARC->delta_y) < 1E-12 && fabs(pSPARC->delta_y - pSPARC->delta_z) < 1E-12) { @@ -1305,51 +1222,14 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) // find effective mesh s.t. it has same spectral width h_eff = sqrt(3.0 / (dx2_inv + dy2_inv + dz2_inv)); } ->>>>>>> variant B - // maximum eigenvalue of -0.5 * Lap (only with periodic boundary conditions) - if(pSPARC->cell_typ == 0) { - #ifdef DEBUG - t1 = MPI_Wtime(); - #endif - pSPARC->MaxEigVal_mhalfLap = pSPARC->D2_stencil_coeffs_x[0] - + pSPARC->D2_stencil_coeffs_y[0] - + pSPARC->D2_stencil_coeffs_z[0]; - double scal_x, scal_y, scal_z; - scal_x = (pSPARC->Nx - pSPARC->Nx % 2) / (double) pSPARC->Nx; - scal_y = (pSPARC->Ny - pSPARC->Ny % 2) / (double) pSPARC->Ny; - scal_z = (pSPARC->Nz - pSPARC->Nz % 2) / (double) pSPARC->Nz; - for (int p = 1; p < FDn + 1; p++) { - pSPARC->MaxEigVal_mhalfLap += 2.0 * (pSPARC->D2_stencil_coeffs_x[p] * cos(M_PI*p*scal_x) - + pSPARC->D2_stencil_coeffs_y[p] * cos(M_PI*p*scal_y) - + pSPARC->D2_stencil_coeffs_z[p] * cos(M_PI*p*scal_z)); - } - pSPARC->MaxEigVal_mhalfLap *= -0.5; - #ifdef DEBUG - t2 = MPI_Wtime(); - if (!rank) printf("Max eigenvalue of -0.5*Lap is %.13f, time taken: %.3f ms\n", - pSPARC->MaxEigVal_mhalfLap, (t2-t1)*1e3); - #endif - } -======= end -<<<<<<< variant A // find Chebyshev polynomial degree based on max eigenvalue (spectral width) if (pSPARC->ChebDegree < 0) { pSPARC->ChebDegree = Mesh2ChebDegree(h_eff); #ifdef DEBUG if (!rank && h_eff < 0.1) { printf("WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enought!\n"); ->>>>>>> variant B - double h_eff = 0.0; - if (fabs(pSPARC->delta_x - pSPARC->delta_y) < 1E-12 && - fabs(pSPARC->delta_y - pSPARC->delta_z) < 1E-12) { - h_eff = pSPARC->delta_x; - } else { - // find effective mesh s.t. it has same spectral width - h_eff = sqrt(3.0 / (dx2_inv + dy2_inv + dz2_inv)); -======= end } -<<<<<<< variant A if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree); #endif } else { @@ -1357,49 +1237,18 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) if (!rank) printf("Chebyshev polynomial degree (provided by user): npl = %d\n",pSPARC->ChebDegree); #endif } ->>>>>>> variant B -======= end -<<<<<<< variant A // default Kerker tolerance if (pSPARC->TOL_PRECOND < 0.0) { // kerker tol not provided by user pSPARC->TOL_PRECOND = (h_eff * h_eff) * 1e-3; } ->>>>>>> variant B - // find Chebyshev polynomial degree based on max eigenvalue (spectral width) - if (pSPARC->ChebDegree < 0) { - pSPARC->ChebDegree = Mesh2ChebDegree(h_eff); - #ifdef DEBUG - if (!rank && h_eff < 0.1) { - printf("#WARNING: for mesh less than 0.1, the default Chebyshev polynomial degree might not be enough!\n"); - } - if (!rank) printf("h_eff = %.2f, npl = %d\n", h_eff,pSPARC->ChebDegree); - #endif - } else { - #ifdef DEBUG - if (!rank) printf("Chebyshev polynomial degree (provided by user): npl = %d\n",pSPARC->ChebDegree); - #endif - } -======= end -<<<<<<< variant A ->>>>>>> variant B - // default Kerker tolerance - if (pSPARC->TOL_PRECOND < 0.0) { // kerker tol not provided by user - pSPARC->TOL_PRECOND = (h_eff * h_eff) * 1e-3; - } - } -======= end // re-calculate k-point grid Calculate_kpoints(pSPARC); // re-calculate local k-points array -<<<<<<< variant A if (pSPARC->Nkpts >= 1 && pSPARC->kptcomm_index != -1) { ->>>>>>> variant B - if (pSPARC->Nkpts >= 1 && pSPARC->kptcomm_index != -1 && (!pSPARC->sqAmbientFlag)) { -======= end Calculate_local_kpoints(pSPARC); } @@ -1413,26 +1262,12 @@ void reinitialize_mesh_NPT(SPARC_OBJ *pSPARC) if (rank == 0) printf("Calculating rb for all atom types took %.3f ms\n",(t2-t1)*1000); #endif -<<<<<<< variant A ->>>>>>> variant B - // Rescale planewaves for spDFT - if (pSPARC->spDFT_Flag == 1) { - spDFT_planewaves(pSPARC); - } -======= end - -<<<<<<< variant A ->>>>>>> variant B -======= end // write reinitialized parameters into output file if (rank == 0) { write_output_reinit_NPT(pSPARC); } -<<<<<<< variant A ->>>>>>> variant B -======= end } @@ -1500,26 +1335,14 @@ void NPT_NP (SPARC_OBJ *pSPARC) { updatePosition(pSPARC); // Reinitialize mesh size and related variables after changing size of cell reinitialize_mesh_NPT(pSPARC); -<<<<<<< variant A ->>>>>>> variant B -======= end // Charge extrapolation (for better rho_guess) elecDensExtrapolation(pSPARC); // Check position of atom near the boundary and apply wraparound in case of PBC, otherwise show error if the atom is too close to the boundary for bounded domain Check_atomlocation(pSPARC); -<<<<<<< variant A - // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem or MLFF - if (pSPARC->mlff_flag == 0){ - Calculate_electronicGroundState(pSPARC); - } else { - MLFF_call(pSPARC); - } ->>>>>>> variant B // Compute DFT energy and forces by solving Kohn-Sham eigenvalue problem Calculate_Properties(pSPARC); //Calculate_electronicGroundState(pSPARC); -======= end pSPARC->elecgs_Count++; #ifdef DEBUG if (!rank) printf("\nend NPT_NP timestep %d\n", pSPARC->MDCount + 1); @@ -1637,12 +1460,6 @@ void updateMomentum_FirstHalf(SPARC_OBJ *pSPARC) { innerControlStress[1] = innerControlStress[0]; innerControlStress[2] = innerControlStress[0]; } -<<<<<<< variant A ->>>>>>> variant B - if (pSPARC->CyclixFlag) { // cyclix on z, unit Ha/Bohr. It is necessary to convert it to Ha/Bohr^3. - innerControlStress[2] /= M_PI * (pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0)) ; - } -======= end // update momentum of barostat variables in a half step for (int i = 0; i < 3; i++){ if (pSPARC->NPTscaleVecs[i] == 1) { @@ -1653,11 +1470,7 @@ void updateMomentum_FirstHalf(SPARC_OBJ *pSPARC) { pSPARC->Pm_NPT_NP[i] -= pSPARC->MD_dt * pSPARC->S_NPT_NP / 2.0 * PmA[i]; #ifdef DEBUG if (rank == 0){ -<<<<<<< variant A // printf("pSPARC->pres is %12.9f, pSPARC->pres_i is %12.9f, pSPARC->volumeCell is %12.9f, pSPARC->G_NPT_NP[%d] is %12.9f\n", pSPARC->pres, pSPARC->pres_i, pSPARC->volumeCell, i, pSPARC->G_NPT_NP[i]); ->>>>>>> variant B - printf("innerControlStress[%d] is %12.9f, pSPARC->volumeCell is %12.9f, pSPARC->G_NPT_NP[%d] is %12.9f\n", i, innerControlStress[i], pSPARC->volumeCell, i, pSPARC->G_NPT_NP[i]); -======= end printf("PmA[%d] is %12.9f\n", i, PmA[i]); printf("pSPARC->Pm_NPT_NP[%d] in 1st half step is %12.9f\n", i, pSPARC->Pm_NPT_NP[i]); } @@ -1721,13 +1534,7 @@ void updateMomentum_SecondHalf(SPARC_OBJ *pSPARC) { innerControlStress[1] = innerControlStress[0]; innerControlStress[2] = innerControlStress[0]; } -<<<<<<< variant A ->>>>>>> variant B - if (pSPARC->CyclixFlag) { // cyclix on z, unit Ha/Bohr. It is necessary to convert it to Ha/Bohr^3. - innerControlStress[2] /= M_PI * (pow(pSPARC->xout,2.0) - pow(pSPARC->xin,2.0)) ; - } -======= end while (judge == 0) { timeIter++; KbaroTmp = 0; @@ -1891,17 +1698,7 @@ void updatePosition(SPARC_OBJ *pSPARC) { pSPARC->range_y = sqrt(pSPARC->G_NPT_NP[1]); double scalez = sqrt(pSPARC->G_NPT_NP[2]) / pSPARC->range_z; pSPARC->range_z = sqrt(pSPARC->G_NPT_NP[2]); -<<<<<<< variant A pSPARC->volumeCell = pSPARC->Jacbdet*pSPARC->range_x*pSPARC->range_y*pSPARC->range_z; ->>>>>>> variant B - if (pSPARC->CyclixFlag) { - pSPARC->twist = pSPARC->twistpercell / pSPARC->range_z; - pSPARC->volumeCell = pSPARC->range_x * ((pSPARC->xin + pSPARC->xout)/2.0) * pSPARC->range_y * pSPARC->range_z; - } - else { - pSPARC->volumeCell = pSPARC->Jacbdet*pSPARC->range_x*pSPARC->range_y*pSPARC->range_z; - } -======= end if (pSPARC->NPTscaleVecs[0] == 1) pSPARC->range_x_velo = G3[0] * Stemp / (2.0*pSPARC->NPT_NP_bmass*pow(pSPARC->volumeCell, 2.0)) / pSPARC->range_x; else @@ -1948,11 +1745,7 @@ void Cart2nonCart(double *gradT, double *carCoord, double *nonCarCoord) { } /* -<<<<<<< variant A * @ brief: function to check if the atoms are too close to the boundary in case of bounded domain or to each other in general ->>>>>>> variant B - @ brief: function to wrap around atom positions that lie outside main domain in case of PBC and check if the atoms are too close to the boundary in case of bounded domain -======= end */ void Check_atomlocation(SPARC_OBJ *pSPARC) { int rank, ityp, i, atm, atm2, count, dir = 0, maxdir = 3, BC; @@ -2000,17 +1793,9 @@ void Check_atomlocation(SPARC_OBJ *pSPARC) { wraparound_dynamics(pSPARC, pSPARC->atom_pos, 1); } -<<<<<<< variant A /* * @ brief: function to wraparound atom positions for PBC */ ->>>>>>> variant B - - - - - -======= end void wraparound_dynamics(SPARC_OBJ *pSPARC, double *coord, int opt) { int rank, atm, dir = 0, maxdir = 3, BC; @@ -2065,13 +1850,9 @@ void wraparound_dynamics(SPARC_OBJ *pSPARC, double *coord, int opt) { } } -<<<<<<< variant A /* * @ brief: function to wraparound velocities in MD and displacement vectors in relaxation for PBC */ ->>>>>>> variant B - -======= end void wraparound_velocity(SPARC_OBJ *pSPARC, double shift, int dir, int loc) { if (dir == 1){ @@ -2100,12 +1881,6 @@ void wraparound_velocity(SPARC_OBJ *pSPARC, double shift, int dir, int loc) { } -<<<<<<< variant A ->>>>>>> variant B - - - -======= end /* @ brief: function to write all relevant DFT quantities generated during MD simulation */ @@ -2355,11 +2130,7 @@ void MD_QOI(SPARC_OBJ *pSPARC, double *avgvel, double *maxvel, double *mindis) { } // Average and minimum distance -<<<<<<< variant A //*avgdis = pow((pSPARC->range_x * pSPARC->range_y * pSPARC->range_z)/pSPARC->n_atom,1/3.0); ->>>>>>> variant B - // *avgdis = pow((pSPARC->range_x * pSPARC->range_y * pSPARC->range_z)/pSPARC->n_atom,1/3.0); -======= end int atm2, ityp2, cc2, pairIndex; cc = 0; for(ityp = 0; ityp < pSPARC->Ntypes; ityp++){ @@ -2490,13 +2261,7 @@ void PrintMD(SPARC_OBJ *pSPARC, int Flag, int print_restart_typ) { fprintf(mdout,":CELL: %.15g %.15g %.15g\n",pSPARC->range_x,pSPARC->range_y,pSPARC->range_z); //(no variable for position of barostat variable) else fprintf(mdout,":LATVEC_SCALE: %.15g %.15g %.15g\n",pSPARC->range_x/pSPARC->initialLatVecLength[0],pSPARC->range_y/pSPARC->initialLatVecLength[1],pSPARC->range_z/pSPARC->initialLatVecLength[2]); -<<<<<<< variant A fprintf(mdout,":TTHRMI(K): %.15g\n", pSPARC->thermos_Ti); ->>>>>>> variant B - if (pSPARC->CyclixFlag) - fprintf(mdout,":TWIST_ANGLE: %.15g\n", pSPARC->twist); - fprintf(mdout,":TTHRMI(K): %.15g\n", pSPARC->thermos_Ti); -======= end fprintf(mdout,":TARGET_PRESSURE: %.15g GPa\n",pSPARC->prtarget * 29421.02648438959); fprintf(mdout,":NPT_NP_ini_Hamiltonian: %.15g\n", pSPARC->init_Hamil_NPT_NP); } @@ -2547,11 +2312,7 @@ void RestartMD(SPARC_OBJ *pSPARC) { l_buff = (2 + 1) * sizeof(int) + (6 * pSPARC->n_atom + (5 + 3*pSPARC->NPT_NHnnos + 9)) * sizeof(double); } else if (strcmpi(pSPARC->MDMeth,"NPT_NP") == 0){ -<<<<<<< variant A l_buff = 2 * sizeof(int) + (6 * pSPARC->n_atom + (5 + 14)) * sizeof(double); ->>>>>>> variant B - l_buff = 2 * sizeof(int) + (6 * pSPARC->n_atom + (5 + 15)) * sizeof(double); -======= end } else { l_buff = 2 * sizeof(int) + (6 * pSPARC->n_atom + 5) * sizeof(double); @@ -2701,11 +2462,6 @@ void RestartMD(SPARC_OBJ *pSPARC) { } else if (strcmpi(str,":TTHRMI(K):") == 0) fscanf(rst_fp,"%lf", &pSPARC->thermos_Ti); -<<<<<<< variant A ->>>>>>> variant B - else if (strcmpi(str,":TWIST_ANGLE:") == 0) - fscanf(rst_fp,"%lf", &pSPARC->twist); -======= end else if (strcmpi(str,":TARGET_PRESSURE:") == 0) fscanf(rst_fp,"%lf", &pSPARC->prtarget); else if (strcmpi(str,":NPT_NP_ini_Hamiltonian:") == 0) @@ -2728,11 +2484,7 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Pack(&pSPARC->xi_nose, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->thermos_Ti, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); } -<<<<<<< variant A else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ ->>>>>>> variant B - else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ -======= end MPI_Pack(&pSPARC->NPT_NHnnos, 1, MPI_INT, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(pSPARC->NPT_NHqmass, pSPARC->NPT_NHnnos, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(pSPARC->vlogs, pSPARC->NPT_NHnnos, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); @@ -2760,10 +2512,6 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Pack(&pSPARC->range_y, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->range_z, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->thermos_Ti, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); -<<<<<<< variant A ->>>>>>> variant B - MPI_Pack(&pSPARC->twist, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); -======= end MPI_Pack(&pSPARC->prtarget, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); MPI_Pack(&pSPARC->init_Hamil_NPT_NP, 1, MPI_DOUBLE, buff, l_buff, &position, MPI_COMM_WORLD); } @@ -2795,11 +2543,7 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Unpack(buff, l_buff, &position, &pSPARC->xi_nose, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->thermos_Ti, 1, MPI_DOUBLE, MPI_COMM_WORLD); } -<<<<<<< variant A else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ ->>>>>>> variant B - else if(strcmpi(pSPARC->MDMeth,"NPT_NH") == 0){ -======= end MPI_Unpack(buff, l_buff, &position, &pSPARC->NPT_NHnnos, 1, MPI_INT, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, pSPARC->NPT_NHqmass, pSPARC->NPT_NHnnos, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, pSPARC->vlogs, pSPARC->NPT_NHnnos, MPI_DOUBLE, MPI_COMM_WORLD); @@ -2828,10 +2572,6 @@ void RestartMD(SPARC_OBJ *pSPARC) { MPI_Unpack(buff, l_buff, &position, &pSPARC->range_y, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->range_z, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->thermos_Ti, 1, MPI_DOUBLE, MPI_COMM_WORLD); -<<<<<<< variant A ->>>>>>> variant B - MPI_Unpack(buff, l_buff, &position, &pSPARC->twist, 1, MPI_DOUBLE, MPI_COMM_WORLD); -======= end MPI_Unpack(buff, l_buff, &position, &pSPARC->prtarget, 1, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Unpack(buff, l_buff, &position, &pSPARC->init_Hamil_NPT_NP, 1, MPI_DOUBLE, MPI_COMM_WORLD); } @@ -2856,7 +2596,4 @@ void Rename_restart(SPARC_OBJ *pSPARC) { if( access(pSPARC->restartC_Filename, F_OK ) != -1 ) rename(pSPARC->restartC_Filename, pSPARC->restartP_Filename); } -<<<<<<< variant A ->>>>>>> variant B -======= end From ee39314d5bd52a0afc126094725f821223908655 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 17:23:50 +0800 Subject: [PATCH 03/14] add use_socket in build --- .conda/build.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.conda/build.sh b/.conda/build.sh index ee9ffeac..3b689811 100644 --- a/.conda/build.sh +++ b/.conda/build.sh @@ -4,7 +4,7 @@ echo Compiling sparc using ${CPU_COUNT} cores cd ./src make clean -make -j ${CPU_COUNT} USE_MKL=0 USE_SCALAPACK=1 +make -j ${CPU_COUNT} USE_MKL=0 USE_SCALAPACK=1 USE_FFTW=1 USE_SOCKET=1 # ls -al lib echo "Installing sparc into $PREFIX/bin" From 3e4e596752c8c9780265953681212d37eb8d3f1e Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 19:07:12 +0800 Subject: [PATCH 04/14] update tests examples --- .github/workflows/build-test.yml | 4 + src/include/isddft.h | 26 +++ src/initialization.c | 114 ++++++++++- tests/Socket/Al_ase_opt/test_all.py | 98 ++++++++++ tests/Socket/Al_geopt/test_all.py | 122 ++++++++++++ tests/Socket/Al_sameimage/test_all.py | 147 ++++++++++++++ tests/Socket/Al_singlepoint_unix/test_all.py | 145 ++++++++++++++ tests/Socket/Al_singlepoints/test_all.py | 141 ++++++++++++++ tests/Socket/Al_volume_change/test_all.py | 150 +++++++++++++++ tests/Socket/Atom_sort_test/test_all.py | 75 ++++++++ tests/Socket/PBC_test/test_all.py | 191 +++++++++++++++++++ tests/Socket/run_all.sh | 12 ++ 12 files changed, 1224 insertions(+), 1 deletion(-) create mode 100644 tests/Socket/Al_ase_opt/test_all.py create mode 100644 tests/Socket/Al_geopt/test_all.py create mode 100644 tests/Socket/Al_sameimage/test_all.py create mode 100644 tests/Socket/Al_singlepoint_unix/test_all.py create mode 100644 tests/Socket/Al_singlepoints/test_all.py create mode 100644 tests/Socket/Al_volume_change/test_all.py create mode 100644 tests/Socket/Atom_sort_test/test_all.py create mode 100644 tests/Socket/PBC_test/test_all.py create mode 100644 tests/Socket/run_all.sh diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 72e2046e..81be871f 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -58,3 +58,7 @@ jobs: unset LD_LIBRARY_PATH python -c "import os, re; match = re.findall(r'Tests failed: (\d+)', open('output').read()); [print('All pass'), os._exit(0)] if int(match[0]) == 0 else [print('Failed'), os._exit(1)]" rm -rf output + - name: Test socket compatibility + run: | + mpirun -n 1 sparc || true + diff --git a/src/include/isddft.h b/src/include/isddft.h index 34d123ea..de0fb238 100644 --- a/src/include/isddft.h +++ b/src/include/isddft.h @@ -1093,6 +1093,20 @@ typedef struct _SPARC_OBJ{ char InDensDCubFilename[L_STRING]; int densfilecount; int readInitDens; // flag for reading inital density + /* Socket interface + Please keep this section as the last block in SPARC_OBJ definition, + add new features before this block. + */ + +#ifdef USE_SOCKET + int SocketFlag; // boolean value to indicate whether to switch on socket support at all + int SocketSCFCount; // scf count in socket mode + char socket_host[L_STRING]; // socket file descriptor + int socket_port; // socket host + int socket_inet; // boolean value to indicate whether to use inet socket + int socket_fd; // socket file descriptor; This should be initialized to -1 + int socket_max_niter; // max number of iterations. Default is 10000 +#endif }SPARC_OBJ; @@ -1404,6 +1418,18 @@ typedef struct _SPARC_INPUT_OBJ{ char InDensDCubFilename[L_STRING]; int densfilecount; int readInitDens; // flag for reading inital density + + /* Socket interface + Please keep the socket interface as the last block in the SPARC_INPUT_OBJ + definition. All new features should be added before this block + */ +#ifdef USE_SOCKET + int SocketFlag; // boolean value to indicate whether to switch on socket support at all + char socket_host[L_STRING]; // socket file descriptor + int socket_port; // socket port + int socket_inet; // boolean value to indicate whether to use inet socket + int socket_max_niter; // max number of iterations +#endif }SPARC_INPUT_OBJ; diff --git a/src/initialization.c b/src/initialization.c index 79ac95d5..c659cf4e 100644 --- a/src/initialization.c +++ b/src/initialization.c @@ -42,6 +42,10 @@ #include "ofdft.h" #include "sparc_mlff_interface.h" +#ifdef USE_SOCKET +#include "driver.h" // socker driver functions +#endif + #define TEMP_TOL 1e-12 #define min(x,y) ((x)<(y)?(x):(y)) @@ -192,6 +196,11 @@ void print_usage() { printf(" -n \n"); printf(" -c \n"); printf(" -a \n"); +#ifdef USE_SOCKET + printf(" -socket \n"); + printf(" can be either : or :UNIX.\n"); + printf(" Note: socket (driver) mode is an experimental feature.\n"); +#endif // USE_SOCKET printf("\n"); printf("EXAMPLE:\n"); printf("\n"); @@ -422,6 +431,11 @@ void Initialize(SPARC_OBJ *pSPARC, int argc, char *argv[]) { } } +#ifdef DEBUG + t2 = MPI_Wtime(); + if (rank == 0) printf("\nrank = %d, Copying data from SPARC_Input into SPARC & set up subcomm took %.3f ms\n",rank,(t2-t1)*1000); +#endif + #ifdef DEBUG t1 = MPI_Wtime(); #endif @@ -492,6 +506,18 @@ void Initialize(SPARC_OBJ *pSPARC, int argc, char *argv[]) { // estimate memory usage pSPARC->memory_usage = estimate_memory(pSPARC); + + + /* NOTE to future developers: Please keep the initialize_Socket as the last + part of the initialization function + */ +#ifdef USE_SOCKET + // Initialize the socket communicator as the last step + if (pSPARC->SocketFlag == 1) + { + initialize_Socket(pSPARC); + } +#endif // write initialized parameters into output file if (rank == 0) { write_output_init(pSPARC); @@ -553,6 +579,29 @@ void check_inputs(SPARC_INPUT_OBJ *pSPARC_Input, int argc, char *argv[], int *ex if (strcmp(argv[i],"-a") == 0) { pSPARC_Input->num_acc_per_node = atoi(argv[i+1]); } +#ifdef USE_SOCKET + if (strcmp(argv[i], "-socket") == 0) + { + const char *socket_str = argv[i + 1]; + pSPARC_Input->SocketFlag = 1; + int ret = split_socket_name(socket_str, + pSPARC_Input->socket_host, + &pSPARC_Input->socket_port, + &pSPARC_Input->socket_inet); + if (ret != 0) + { + printf("Error: invalid socket name %s\n", socket_str); + exit(EXIT_FAILURE); + } +#ifdef DEBUG + printf("Socket host = %s\n", pSPARC_Input->socket_host); + printf("Socket port is %d\n", pSPARC_Input->socket_port); + printf("Socket inet flag is %d\n", pSPARC_Input->socket_inet); + printf("Socket mode is %d\n", pSPARC_Input->SocketFlag); +#endif // DEBUG + } + +#endif // USE_SOCKET } if (name_flag != 'Y') { @@ -857,6 +906,31 @@ void set_defaults(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) { // read in initial density pSPARC_Input->readInitDens = 0; + + /* Default socket options + Note to future developers: please keep the USE_SOCKET macro + as the last part of the initialization function + */ +#ifdef USE_SOCKET + // Defaults for pSPARC_Input, if not initialized by the cmdline + if (pSPARC_Input->SocketFlag != 1) + { + pSPARC_Input->SocketFlag = 0; // socket off + strncpy(pSPARC_Input->socket_host, "localhost", sizeof(pSPARC_Input->socket_host)); + pSPARC_Input->socket_port = -1; // socket port + pSPARC_Input->socket_inet = 0; // 0: -> default unix socket, 1: -> inet socket + } + pSPARC_Input->socket_max_niter = 10000; // Set to a very large number +#endif + + + /* Default spDFT option */ + pSPARC_Input->spDFT_Flag = 0; // flag for using spDFT method + pSPARC_Input->PrintspDFTFlag = 0; // flag for printing spDFT data into .spDFT file + pSPARC_Input->spDFT_isesplit_const = 0; // Is splitting energy to be kept constant + pSPARC_Input->spDFT_tau_s = 0.01; // default smearing (in Ha) for spectral-partitioning + pSPARC_Input->spDFT_tol_occ = 1e-12; // default maximum occupation of the highest-energy planewave + } @@ -1461,6 +1535,25 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) { strncpy(pSPARC->InDensTCubFilename, pSPARC_Input->InDensTCubFilename,sizeof(pSPARC->InDensTCubFilename)); strncpy(pSPARC->InDensUCubFilename, pSPARC_Input->InDensUCubFilename,sizeof(pSPARC->InDensUCubFilename)); strncpy(pSPARC->InDensDCubFilename, pSPARC_Input->InDensDCubFilename,sizeof(pSPARC->InDensDCubFilename)); + + /* Socket interface section + TODO: should we move the socket to a later section? + */ +#ifdef USE_SOCKET + // Copy the pSPARC_Input socket information to pSPARC. + // Since only rank 0 handles the communication, we don't need + // to explicitly bcast except for SocketFlag + pSPARC->SocketFlag = pSPARC_Input->SocketFlag; + pSPARC->socket_inet = pSPARC_Input->socket_inet; + strncpy(pSPARC->socket_host, pSPARC_Input->socket_host, sizeof(pSPARC->socket_host)); + pSPARC->socket_port = pSPARC_Input->socket_port; + pSPARC->socket_max_niter = pSPARC_Input->socket_max_niter; + // pSPARC->SocketSCFCount needs to be initialized here for non-socket calculations + pSPARC->SocketSCFCount = 0; + // Only SocketFlag and socket_max_ninter information is meaningful at all ranks + MPI_Bcast(&pSPARC->SocketFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&pSPARC->socket_max_niter, 1, MPI_INT, 0, MPI_COMM_WORLD); +#endif if (pSPARC->BandStructFlag == 1) { if (pSPARC->densfilecount != 3 && pSPARC->spin_typ == 1) { @@ -3564,7 +3657,7 @@ void write_output_init(SPARC_OBJ *pSPARC) { } fprintf(output_fp,"***************************************************************************\n"); - fprintf(output_fp,"* SPARC (version Oct 30, 2024) *\n"); + fprintf(output_fp,"* SPARC (version Nov 18, 2024) *\n"); fprintf(output_fp,"* Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech *\n"); fprintf(output_fp,"* Distributed under GNU General Public License 3 (GPL) *\n"); fprintf(output_fp,"* Start time: %s *\n",c_time_str); @@ -3972,6 +4065,20 @@ void write_output_init(SPARC_OBJ *pSPARC) { fprintf(output_fp,"OUTPUT_FILE: %s\n",pSPARC->filename_out); + +#ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) + { + fprintf(output_fp, "***************************************************************************\n"); + fprintf(output_fp, " Socket Mode \n"); + fprintf(output_fp, "***************************************************************************\n"); + fprintf(output_fp, "SOCKET_HOST: %s\n", pSPARC->socket_host); + fprintf(output_fp, "SOCKET_PORT: %d\n", pSPARC->socket_port); + fprintf(output_fp, "SOCKET_INET: %d\n", pSPARC->socket_inet); + fprintf(output_fp, "SOCKET_MAX_NITER: %d\n", pSPARC->socket_max_niter); + } +#endif // USE_SOCKET + fprintf(output_fp,"***************************************************************************\n"); fprintf(output_fp," Cell \n"); fprintf(output_fp,"***************************************************************************\n"); @@ -4093,6 +4200,7 @@ void write_output_init(SPARC_OBJ *pSPARC) { fclose(output_fp); // write .static file +#ifndef USE_SOCKET if ((pSPARC->PrintAtomPosFlag == 1 || pSPARC->PrintForceFlag == 1) && pSPARC->MDFlag == 0 && pSPARC->RelaxFlag == 0) { FILE *static_fp = fopen(pSPARC->StaticFilename,"w"); if (static_fp == NULL) { @@ -4119,6 +4227,10 @@ void write_output_init(SPARC_OBJ *pSPARC) { } fclose(static_fp); } +#else + // Use the print method in socket driver to print the static file + socket_static_print_atom_pos(pSPARC); +#endif } diff --git a/tests/Socket/Al_ase_opt/test_all.py b/tests/Socket/Al_ase_opt/test_all.py new file mode 100644 index 00000000..eb6c92b0 --- /dev/null +++ b/tests/Socket/Al_ase_opt/test_all.py @@ -0,0 +1,98 @@ +"""Testing single point calculations between pure SPARC and socket mode +""" + +import os +import shutil +import time +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write +from ase.optimize.lbfgs import LBFGS +from sparc import SPARC + +os.environ["SPARC_PSP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.28, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + "RELAX_FLAG": 1, + "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + + +al = bulk("Al", cubic=True) +al.rattle(0.2, seed=42) + + +def sparc_singlepoint(): + shutil.rmtree("sparc_geopt", ignore_errors=True) + try: + os.remove("sparc_geopt.extxyz") + except Exception: + pass + atoms = al.copy() + calc = SPARC( + directory=f"sparc_geopt", + command="mpirun -n 2 ../../../../lib/sparc", + **sparc_params, + ) + atoms.calc = calc + e = atoms.get_potential_energy() + out_images = read("sparc_geopt", format="sparc", index=":") + write("sparc_geopt.extxyz", out_images) + return out_images + + +def sparc_socket(): + atoms = al.copy() + + inputs = Path("./sparc_geopt") + copy_to = Path("./socket_test") + shutil.rmtree(copy_to, ignore_errors=True) + os.makedirs(copy_to, exist_ok=True) + shutil.copy(inputs / "SPARC.ion", copy_to) + shutil.copy(inputs / "SPARC.inpt", copy_to) + shutil.copy(inputs / "13_Al_3_1.9_1.9_pbe_n_v1.0.psp8", copy_to) + + with SocketIOCalculator(port=12345) as calc: + time.sleep( + 1.0 + ) # In some emulated systems the SocketIOCalculator may delay binding + p_ = Popen( + "mpirun -n 2 ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log", + shell=True, + cwd=copy_to, + ) + atoms.calc = calc + opt = LBFGS(atoms, trajectory="sparc-socket.traj") + opt.run(fmax=0.02) + p_.kill() + return atoms.copy() + + +def main(): + images_sparc = sparc_singlepoint() + final_socket = sparc_socket() + final_sparc = images_sparc[-1] + final_socket.wrap() + final_sparc.wrap() + positions_change = final_socket.positions - final_sparc.positions + max_change = np.linalg.norm(positions_change) + print("Max shift: ", max_change) + assert max_change < 0.02 + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/Al_geopt/test_all.py b/tests/Socket/Al_geopt/test_all.py new file mode 100644 index 00000000..0b3f3378 --- /dev/null +++ b/tests/Socket/Al_geopt/test_all.py @@ -0,0 +1,122 @@ +"""Testing single point calculations between pure SPARC and socket mode +""" + +import os +import shutil +import time +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write +from sparc import SPARC + +os.environ["SPARC_PP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.28, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + "RELAX_FLAG": 1, + "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + + +al = bulk("Al", cubic=True) +al.rattle(0.2, seed=42) + + +def sparc_singlepoint(): + shutil.rmtree("sparc_geopt", ignore_errors=True) + try: + os.remove("sparc_geopt.extxyz") + except Exception: + pass + atoms = al.copy() + calc = SPARC( + directory=f"sparc_geopt", + command="mpirun -n 2 --oversubscribe ../../../../lib/sparc", + **sparc_params, + ) + atoms.calc = calc + e = atoms.get_potential_energy() + out_images = read("sparc_geopt", format="sparc", index=":") + write("sparc_geopt.extxyz", out_images) + return out_images + + +def sparc_socket(): + imgfile = Path("sparc_geopt.extxyz") + images = read(imgfile, ":") + + inputs = Path("./sparc_geopt") + copy_to = Path("./socket_test") + shutil.rmtree(copy_to, ignore_errors=True) + os.makedirs(copy_to, exist_ok=True) + shutil.copy(inputs / "SPARC.ion", copy_to) + shutil.copy(inputs / "SPARC.inpt", copy_to) + shutil.copy(inputs / "13_Al_3_1.9_1.9_pbe_n_v1.0.psp8", copy_to) + + out_images = [] + with SocketIOCalculator(port=12345) as calc: + time.sleep(1.0) + p_ = Popen( + "mpirun -n 2 --oversubscribe ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log", + shell=True, + cwd=copy_to, + ) + for i, atoms in enumerate(images): + atoms.calc = calc + e = atoms.get_potential_energy() + f = atoms.get_forces() + s = atoms.get_stress() + atoms_cp = atoms.copy() + atoms_cp.calc = SinglePointCalculator( + atoms_cp, energy=e, forces=f, stress=s + ) + out_images.append(atoms_cp) + p_.kill() + return out_images + + +def main(): + images_sparc = sparc_singlepoint() + images_socket = sparc_socket() + ediff = np.array( + [ + img2.get_potential_energy() - img1.get_potential_energy() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + fdiff = np.array( + [ + img2.get_forces() - img1.get_forces() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + sdiff = np.array( + [ + img2.get_stress() - img1.get_stress() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + + print("Ediff, ", np.max(np.abs(ediff))) + print("Fdiff, ", np.max(np.abs(fdiff))) + print("Sdiff, ", np.max(np.abs(sdiff))) + assert np.max(np.abs(ediff)) < 1.0e-3 + assert np.max(np.abs(fdiff)) < 0.01 + assert np.max(np.abs(sdiff)) < 0.01 + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/Al_sameimage/test_all.py b/tests/Socket/Al_sameimage/test_all.py new file mode 100644 index 00000000..4c9c781e --- /dev/null +++ b/tests/Socket/Al_sameimage/test_all.py @@ -0,0 +1,147 @@ +"""Testing single point calculations between pure SPARC and socket mode +""" + +import os +import shutil +import time +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write +from sparc import SPARC + +os.environ["SPARC_PSP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.28, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + "RELAX_FLAG": 0, + "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + + +def make_images(): + atoms = bulk("Al", cubic=True) + + imgfile = Path("al_images.extxyz") + + if not imgfile.is_file(): + images = [] + for i in range(5): + at = atoms.copy() + at.rattle(0.1, seed=i) + images.append(at) + write(imgfile, images) + + images = read(imgfile, ":") + return images + + +def sparc_singlepoint(): + for d in Path(".").glob("sp_image*"): + shutil.rmtree(d, ignore_errors=True) + images = make_images()[:1] + out_images = [] + for i, img in enumerate(images): + # relateive to the path + calc = SPARC( + directory=f"sp_image{i:02d}", + command="mpirun -n 2 --oversubscribe ../../../../lib/sparc", + **sparc_params, + ) + img.calc = calc + e = img.get_potential_energy() + f = img.get_forces() + s = img.get_stress() + img_cp = img.copy() + img_cp.calc = SinglePointCalculator(img_cp, energy=e, forces=f, stress=s) + out_images.append(img_cp) + print(f"image {i}, energy {e}") + return out_images + + +def sparc_socket(): + images = make_images()[:1] + # Copy inputs from single point --> socket + inputs = Path("./sp_image00") + copy_to = Path("./socket_test") + shutil.rmtree(copy_to, ignore_errors=True) + os.makedirs(copy_to, exist_ok=True) + shutil.copy(inputs / "SPARC.ion", copy_to) + shutil.copy(inputs / "SPARC.inpt", copy_to) + shutil.copy(inputs / "13_Al_3_1.9_1.9_pbe_n_v1.0.psp8", copy_to) + + atoms = images[0] + out_images = [] + with SocketIOCalculator(port=12345) as calc: + time.sleep(1) + p_ = Popen( + "mpirun -n 2 --oversubscribe ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log 2>&1", + shell=True, + cwd=copy_to, + ) + for i in range(5): + atoms.calc = calc + # Force clean results + calc.results.clear() + e = atoms.get_potential_energy() + f = atoms.get_forces() + s = atoms.get_stress() + e_old = read("sp_image00", format="sparc", index=-1).get_potential_energy() + # forces = atoms.get_forces() + # stress = atoms.get_stress() + print("Cycle: ", i) + print("Energy: ", e) + print("Energy SPARC SP: ", e_old) + atoms_cp = atoms.copy() + atoms_cp.calc = SinglePointCalculator( + atoms_cp, energy=e, forces=f, stress=s + ) + out_images.append(atoms_cp) + p_.kill() + return out_images + + +def main(): + images_sparc = sparc_singlepoint() + images_socket = sparc_socket() + ediff = np.array( + [ + img2.get_potential_energy() - img1.get_potential_energy() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + fdiff = np.array( + [ + img2.get_forces() - img1.get_forces() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + sdiff = np.array( + [ + img2.get_stress() - img1.get_stress() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + + print("Ediff, ", np.max(np.abs(ediff))) + print("Fdiff, ", np.max(np.abs(fdiff))) + print("Sdiff, ", np.max(np.abs(sdiff))) + assert np.max(np.abs(ediff)) < 1.0e-8 + assert np.max(np.abs(fdiff)) < 1.0e-8 + assert np.max(np.abs(sdiff)) < 1.0e-8 + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/Al_singlepoint_unix/test_all.py b/tests/Socket/Al_singlepoint_unix/test_all.py new file mode 100644 index 00000000..6db74ea0 --- /dev/null +++ b/tests/Socket/Al_singlepoint_unix/test_all.py @@ -0,0 +1,145 @@ +"""Testing single point calculations between pure SPARC and socket mode + specifically use the Unix socket +""" + +import os +import shutil +import time +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write +from sparc import SPARC + +os.environ["SPARC_PP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.28, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + "RELAX_FLAG": 0, + "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + + +def make_images(): + atoms = bulk("Al", cubic=True) + + imgfile = Path("al_images.extxyz") + + if not imgfile.is_file(): + images = [] + for i in range(5): + at = atoms.copy() + at.rattle(0.1, seed=i) + images.append(at) + write(imgfile, images) + + images = read(imgfile, ":") + return images + + +def sparc_singlepoint(): + for d in Path(".").glob("sp_image*"): + shutil.rmtree(d, ignore_errors=True) + images = make_images() + out_images = [] + for i, img in enumerate(images): + # relateive to the path + calc = SPARC( + directory=f"sp_image{i:02d}", + command="mpirun -n 2 --oversubscribe ../../../../lib/sparc", + **sparc_params, + ) + img.calc = calc + e = img.get_potential_energy() + f = img.get_forces() + s = img.get_stress() + img_cp = img.copy() + img_cp.calc = SinglePointCalculator(img_cp, energy=e, forces=f, stress=s) + out_images.append(img_cp) + print(f"image {i}, energy {e}") + return out_images + + +def sparc_socket(): + images = make_images() + # Copy inputs from single point --> socket + inputs = Path("./sp_image00") + copy_to = Path("./socket_test") + shutil.rmtree(copy_to, ignore_errors=True) + os.makedirs(copy_to, exist_ok=True) + shutil.copy(inputs / "SPARC.ion", copy_to) + shutil.copy(inputs / "SPARC.inpt", copy_to) + shutil.copy(inputs / "13_Al_3_1.9_1.9_pbe_n_v1.0.psp8", copy_to) + + out_images = [] + with SocketIOCalculator(unixsocket="sparc") as calc: + time.sleep(1.0) + p_ = Popen( + "mpirun -n 2 --oversubscribe ../../../../lib/sparc -socket /tmp/ipi_sparc:unix -name SPARC > sparc.log 2>&1", + shell=True, + cwd=copy_to, + ) + for i, atoms in enumerate(images): + atoms.calc = calc + e = atoms.get_potential_energy() + f = atoms.get_forces() + s = atoms.get_stress() + e_old = read( + f"sp_image{i:02d}", format="sparc", index=-1 + ).get_potential_energy() + print("Cycle: ", i) + print("Energy: ", e) + print("Energy SPARC SP: ", e_old) + atoms_cp = atoms.copy() + atoms_cp.calc = SinglePointCalculator( + atoms_cp, energy=e, forces=f, stress=s + ) + out_images.append(atoms_cp) + p_.kill() + return out_images + + +def main(): + images_sparc = sparc_singlepoint() + images_socket = sparc_socket() + ediff = np.array( + [ + img2.get_potential_energy() - img1.get_potential_energy() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + fdiff = np.array( + [ + img2.get_forces() - img1.get_forces() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + sdiff = np.array( + [ + img2.get_stress() - img1.get_stress() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + + print("Ediff, ", np.max(np.abs(ediff))) + print("Fdiff, ", np.max(np.abs(fdiff))) + print("Sdiff, ", np.max(np.abs(sdiff))) + assert np.max(np.abs(ediff)) < 1.0e-4 + assert np.max(np.abs(fdiff)) < 0.01 + assert np.max(np.abs(sdiff)) < 1.0e-3 + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/Al_singlepoints/test_all.py b/tests/Socket/Al_singlepoints/test_all.py new file mode 100644 index 00000000..1f68385d --- /dev/null +++ b/tests/Socket/Al_singlepoints/test_all.py @@ -0,0 +1,141 @@ +"""Testing single point calculations between pure SPARC and socket mode +""" + +import os +import shutil +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write +from sparc import SPARC + +os.environ["SPARC_PP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.28, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + "RELAX_FLAG": 0, + "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + + +def make_images(): + atoms = bulk("Al", cubic=True) + + imgfile = Path("al_images.extxyz") + + if not imgfile.is_file(): + images = [] + for i in range(5): + at = atoms.copy() + at.rattle(0.1, seed=i) + images.append(at) + write(imgfile, images) + + images = read(imgfile, ":") + return images + + +def sparc_singlepoint(): + for d in Path(".").glob("sp_image*"): + shutil.rmtree(d, ignore_errors=True) + images = make_images() + out_images = [] + for i, img in enumerate(images): + # relateive to the path + calc = SPARC( + directory=f"sp_image{i:02d}", + command="mpirun -n 2 --oversubscribe ../../../../lib/sparc", + **sparc_params, + ) + img.calc = calc + e = img.get_potential_energy() + f = img.get_forces() + s = img.get_stress() + img_cp = img.copy() + img_cp.calc = SinglePointCalculator(img_cp, energy=e, forces=f, stress=s) + out_images.append(img_cp) + print(f"image {i}, energy {e}") + return out_images + + +def sparc_socket(): + images = make_images() + # Copy inputs from single point --> socket + inputs = Path("./sp_image00") + copy_to = Path("./socket_test") + shutil.rmtree(copy_to, ignore_errors=True) + os.makedirs(copy_to, exist_ok=True) + shutil.copy(inputs / "SPARC.ion", copy_to) + shutil.copy(inputs / "SPARC.inpt", copy_to) + shutil.copy(inputs / "13_Al_3_1.9_1.9_pbe_n_v1.0.psp8", copy_to) + + calc = SocketIOCalculator(port=12345) + p_ = Popen( + "mpirun -n 2 --oversubscribe ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log 2>&1", + shell=True, + cwd=copy_to, + ) + out_images = [] + with calc: + for i, atoms in enumerate(images): + atoms.calc = calc + e = atoms.get_potential_energy() + f = atoms.get_forces() + s = atoms.get_stress() + e_old = read( + f"sp_image{i:02d}", format="sparc", index=-1 + ).get_potential_energy() + # forces = atoms.get_forces() + # stress = atoms.get_stress() + print("Cycle: ", i) + print("Energy: ", e) + print("Energy SPARC SP: ", e_old) + atoms_cp = atoms.copy() + atoms_cp.calc = SinglePointCalculator( + atoms_cp, energy=e, forces=f, stress=s + ) + out_images.append(atoms_cp) + return out_images + + +def main(): + images_sparc = sparc_singlepoint() + images_socket = sparc_socket() + ediff = np.array( + [ + img2.get_potential_energy() - img1.get_potential_energy() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + fdiff = np.array( + [ + img2.get_forces() - img1.get_forces() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + sdiff = np.array( + [ + img2.get_stress() - img1.get_stress() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + + print("Ediff, ", np.max(np.abs(ediff))) + print("Fdiff, ", np.max(np.abs(fdiff))) + print("Sdiff, ", np.max(np.abs(sdiff))) + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/Al_volume_change/test_all.py b/tests/Socket/Al_volume_change/test_all.py new file mode 100644 index 00000000..128b8974 --- /dev/null +++ b/tests/Socket/Al_volume_change/test_all.py @@ -0,0 +1,150 @@ +"""Testing single point calculations with volume and position changes +between pure SPARC and socket mode +""" + +import os +import shutil +import time +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write +from sparc import SPARC + +os.environ["SPARC_PSP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.22, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + "RELAX_FLAG": 0, + "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + +rats = [1.0, 1.01, 1.02, 1.03, 1.04, 1.0, 0.99, 0.98, 0.97, 0.96] + +def make_images(): + atoms = bulk("Al", cubic=True) + + images = [] + for i in range(len(rats)): + at = atoms.copy() + # at.rattle(0.1, seed=i) + cell_origin = at.cell.copy() + rat = rats[i] + cell = cell_origin * [rat, rat, rat] + at.set_cell(cell, scale_atoms=True) + images.append(at) + + return images + + +def sparc_singlepoint(): + for d in Path(".").glob("sp_image*"): + shutil.rmtree(d, ignore_errors=True) + images = make_images() + out_images = [] + for i, img in enumerate(images): + # relateive to the path + rat = rats[i] + params = sparc_params.copy() + params["h"] *= rat + calc = SPARC( + directory=f"sp_image{i:02d}", + command="mpirun -n 2 --oversubscribe ../../../../lib/sparc", + **params, + ) + img.calc = calc + e = img.get_potential_energy() + f = img.get_forces() + s = img.get_stress() + img_cp = img.copy() + img_cp.calc = SinglePointCalculator(img_cp, energy=e, forces=f, stress=s) + out_images.append(img_cp) + print(f"image {i}, energy {e}") + return out_images + + +def sparc_socket(): + images = make_images() + # Copy inputs from single point --> socket + inputs = Path("./sp_image00") + copy_to = Path("./socket_test") + shutil.rmtree(copy_to, ignore_errors=True) + os.makedirs(copy_to, exist_ok=True) + shutil.copy(inputs / "SPARC.ion", copy_to) + shutil.copy(inputs / "SPARC.inpt", copy_to) + shutil.copy(inputs / "13_Al_3_1.9_1.9_pbe_n_v1.0.psp8", copy_to) + + out_images = [] + with SocketIOCalculator(port=12345) as calc: + time.sleep(1.0) + p_ = Popen( + "mpirun -n 2 --oversubscribe ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log 2>&1", + shell=True, + cwd=copy_to, + ) + for i, atoms in enumerate(images): + atoms.calc = calc + e = atoms.get_potential_energy() + f = atoms.get_forces() + s = atoms.get_stress() + e_old = read( + f"sp_image{i:02d}", format="sparc", index=-1 + ).get_potential_energy() + forces = atoms.get_forces() + stress = atoms.get_stress() + print("Cycle: ", i) + print("Energy: ", e) + print("Energy SPARC SP: ", e_old) + print("Ediff: ", np.abs(e - e_old)) + atoms_cp = atoms.copy() + atoms_cp.calc = SinglePointCalculator( + atoms_cp, energy=e, forces=f, stress=s + ) + out_images.append(atoms_cp) + p_.kill() + return out_images + + +def main(): + images_sparc = sparc_singlepoint() + images_socket = sparc_socket() + ediff = np.array( + [ + img2.get_potential_energy() - img1.get_potential_energy() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + fdiff = np.array( + [ + img2.get_forces() - img1.get_forces() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + sdiff = np.array( + [ + img2.get_stress() - img1.get_stress() + for img1, img2 in zip(images_sparc, images_socket) + ] + ) + + print("Ediff, ", np.max(np.abs(ediff))) + print("Fdiff, ", np.max(np.abs(fdiff))) + print("Sdiff, ", np.max(np.abs(sdiff))) + assert np.max(np.abs(ediff)) < 1.e-3 + assert np.max(np.abs(fdiff)) < 5.e-3 + assert np.max(np.abs(ediff)) < 1.e-3 + +if __name__ == "__main__": + main() diff --git a/tests/Socket/Atom_sort_test/test_all.py b/tests/Socket/Atom_sort_test/test_all.py new file mode 100644 index 00000000..263ceee5 --- /dev/null +++ b/tests/Socket/Atom_sort_test/test_all.py @@ -0,0 +1,75 @@ +"""Testing single point calculations between pure SPARC and socket mode +""" + +import os +import shutil +import time +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import molecule +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write + +# from ase.optimize.lbfgs import LBFGS +from ase.optimize.bfgs import BFGS +from sparc import SPARC + +os.environ["SPARC_PSP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.16, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + # "RELAX_FLAG": 1, + # "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + +# The atoms are sorted +mol = molecule("H2O", vacuum=6, pbc=False)[[1, 2, 0]] + + +def simple_sparc(): + atoms = mol.copy() + calc = SPARC(directory="pbc_sp", **sparc_params) + atoms.calc = calc + opt = BFGS(atoms, trajectory="sparc-sp.traj") + opt.run(fmax=0.02) + return + + +def sparc_socket(): + atoms = mol.copy() + calc = SPARC(directory="h2o_test", **sparc_params) + calc.write_input(atoms) + + with SocketIOCalculator(port=12345) as calc: + time.sleep( + 1.0 + ) # In some emulated systems the SocketIOCalculator may delay binding + p_ = Popen( + "mpirun -n 2 ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log", + shell=True, + cwd="h2o_test", + ) + atoms.calc = calc + opt = BFGS(atoms, trajectory="sparc-socket.traj") + opt.run(fmax=0.03) + assert opt.get_number_of_steps() <= 10 # For the correct mesh spacing, opt should be very close + return + + +def main(): + sparc_socket() + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/PBC_test/test_all.py b/tests/Socket/PBC_test/test_all.py new file mode 100644 index 00000000..7bc62efe --- /dev/null +++ b/tests/Socket/PBC_test/test_all.py @@ -0,0 +1,191 @@ +"""Testing single point calculations between pure SPARC and socket mode +""" + +import os +import shutil +import signal +import time +from contextlib import contextmanager +from pathlib import Path +from subprocess import PIPE, Popen + +import ase +import numpy as np +from ase.build import bulk, molecule +from ase.calculators.singlepoint import SinglePointCalculator +from ase.calculators.socketio import SocketIOCalculator +from ase.io import read, write + +# from ase.optimize.lbfgs import LBFGS +from ase.optimize.bfgs import BFGS +from sparc import SPARC + +os.environ["SPARC_PSP_PATH"] = "../../../psps/" + +sparc_params = { + "h": 0.28, + "PRECOND_KERKER_THRESH": 0, + "ELEC_TEMP_TYPE": "Fermi-Dirac", + "ELEC_TEMP": 300, + "MIXING_PARAMETER": 1.0, + "TOL_SCF": 1e-3, + # "RELAX_FLAG": 1, + # "CALC_STRESS": 1, + "PRINT_ATOMS": 1, + "PRINT_FORCES": 1, +} + + +class TimeoutException(Exception): + """Simple class for timeout""" + + pass + + +@contextmanager +def time_limit(seconds): + """Usage: + try: + with time_limit(60): + do_something() + except TimeoutException: + raise + """ + + def signal_handler(signum, frame): + raise TimeoutException("Timed out closing sparc process.") + + signal.signal(signal.SIGALRM, signal_handler) + signal.alarm(seconds) + try: + yield + finally: + signal.alarm(0) + + +mol_orign = molecule("H2", vacuum=3, pbc=True) + + +def sparc_socket(mol, calc_stress=False): + atoms = mol.copy() + params = sparc_params.copy() + if calc_stress: + params["CALC_STRESS"] = 1 + else: + params["CALC_STRESS"] = 0 + calc_origin = SPARC(directory="pbc_test", **params) + calc_origin.write_input(atoms) + + with SocketIOCalculator(port=12345) as calc: + time.sleep( + 1.0 + ) # In some emulated systems the SocketIOCalculator may delay binding + p_ = Popen( + "mpirun -n 2 ../../../../lib/sparc -socket :12345 -name SPARC > sparc.log", + shell=True, + cwd="pbc_test", + ) + atoms.calc = calc + calc.calculate(atoms) + results = calc.results + return results + + +def main(): + mol = mol_orign.copy() + mol.pbc = True + # CASE 1: pbc = True, there is a stress (although meaningless) + res1 = sparc_socket(mol, calc_stress=True) + print(res1) + assert "stress" in res1.keys() + for i in (0, 1, 2): + assert abs(res1["stress"][i]) > 1.0e-6 # Should see a real stress component + assert abs(res1["stress"][i + 3]) < 1.0e-6 # Should see a real stress component + + # CASE 2: pbc = True but no stress calculated + mol = mol_orign.copy() + mol.pbc = True + res2 = sparc_socket(mol, calc_stress=False) + print(res2) + assert "stress" in res2.keys() + assert np.isclose(res2["stress"], 0).all() + + # CASE 2-1: pbc = False with CALC_STRESS=1 --> cannot continue SPARC calculation + # mol = mol_orign.copy() + # mol.pbc = False + + # CASE 2-2: pbc = False with CALC_STRESS=0. SocketServer will not return stress + # because pbc = False in all directions + mol = mol_orign.copy() + mol.pbc = False + res2_2 = sparc_socket(mol, calc_stress=False) + print(res2_2) + assert "stress" not in res2_2.keys() + + # CASE 3: pbc=True, CALC_STRESS=1 + al = bulk("Al", cubic=True) + res3 = sparc_socket(al, calc_stress=True) + assert "stress" in res3.keys() + for i in (0, 1, 2): + assert abs(res3["stress"][i]) > 1.0e-6 # Should see a real stress component + assert abs(res3["stress"][i + 3]) < 1.0e-6 # Should see a real stress component + + # CASE 4: pbc=[T, T, F], CALC_STRESS=1, and all other 2D BCs + al = bulk("Al", cubic=True) + al.pbc = [True, True, False] + res4 = sparc_socket(al, calc_stress=True) + print(res4) + assert "stress" in res4.keys() + for i in (0, 1): + assert abs(res4["stress"][i]) > 1.0e-6 # Should see a real stress component + assert np.isclose(res4["stress"][2:], 0, atol=1.0e-6).all() + + al = bulk("Al", cubic=True) + al.pbc = [True, False, True] + res4_1 = sparc_socket(al, calc_stress=True) + print(res4_1) + assert "stress" in res4_1.keys() + for i in (0, 2): + assert abs(res4_1["stress"][i]) > 1.0e-6 # Should see a real stress component + assert np.isclose(res4_1["stress"][[1, 3, 4, 5]], 0, atol=1.0e-6).all() + + al = bulk("Al", cubic=True) + al.pbc = [False, True, True] + res4_2 = sparc_socket(al, calc_stress=True) + print(res4_2) + assert "stress" in res4_2.keys() + for i in (1, 2): + assert abs(res4_2["stress"][i]) > 1.0e-6 # Should see a real stress component + assert np.isclose(res4_2["stress"][[0, 3, 4, 5]], 0, atol=1.0e-6).all() + + # CASE 5: pbc=[T, F, F], CALC_STRESS=1, + al = bulk("Al", cubic=True) + al.pbc = [True, False, False] + res5 = sparc_socket(al, calc_stress=True) + print(res5) + assert "stress" in res5.keys() + i = 0 + assert abs(res5["stress"][i]) > 1.0e-6 # Should see a real stress component + assert np.isclose(res5["stress"][1:], 0, atol=1.0e-6).all() + + al = bulk("Al", cubic=True) + al.pbc = [False, True, False] + res5_1 = sparc_socket(al, calc_stress=True) + print(res5_1) + assert "stress" in res5_1.keys() + i = 1 + assert abs(res5_1["stress"][i]) > 1.0e-6 # Should see a real stress component + assert np.isclose(res5_1["stress"][[0, 2, 3, 4, 5]], 0, atol=1.0e-6).all() + + al = bulk("Al", cubic=True) + al.pbc = [False, False, True] + res5_2 = sparc_socket(al, calc_stress=True) + print(res5_2) + assert "stress" in res5_2.keys() + i = 2 + assert abs(res5_2["stress"][i]) > 1.0e-6 # Should see a real stress component + assert np.isclose(res5_2["stress"][[0, 1, 3, 4, 5]], 0, atol=1.0e-6).all() + + +if __name__ == "__main__": + main() diff --git a/tests/Socket/run_all.sh b/tests/Socket/run_all.sh new file mode 100644 index 00000000..fc45274d --- /dev/null +++ b/tests/Socket/run_all.sh @@ -0,0 +1,12 @@ +#!/bin/bash +PWD=$(pwd) +for dir in $PWD/*/ +do + dir=${dir%*/} + + if [[ -f "$dir/test_all.py" ]]; then + echo "Running test_all.py in $dir" + # Change into the directory and run the script + (cd "$dir" && python test_all.py) + fi +done From a94d0a1b307f59b22d234c354961a41945e67205 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 19:23:05 +0800 Subject: [PATCH 05/14] remove spDFT from initialization.c --- src/initialization.c | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/initialization.c b/src/initialization.c index c659cf4e..d56b49cb 100644 --- a/src/initialization.c +++ b/src/initialization.c @@ -909,7 +909,7 @@ void set_defaults(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) { /* Default socket options Note to future developers: please keep the USE_SOCKET macro - as the last part of the initialization function + as the LAST PART of the initialization function!! */ #ifdef USE_SOCKET // Defaults for pSPARC_Input, if not initialized by the cmdline @@ -923,14 +923,6 @@ void set_defaults(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) { pSPARC_Input->socket_max_niter = 10000; // Set to a very large number #endif - - /* Default spDFT option */ - pSPARC_Input->spDFT_Flag = 0; // flag for using spDFT method - pSPARC_Input->PrintspDFTFlag = 0; // flag for printing spDFT data into .spDFT file - pSPARC_Input->spDFT_isesplit_const = 0; // Is splitting energy to be kept constant - pSPARC_Input->spDFT_tau_s = 0.01; // default smearing (in Ha) for spectral-partitioning - pSPARC_Input->spDFT_tol_occ = 1e-12; // default maximum occupation of the highest-energy planewave - } From 9baf1e53830434728b2b1bfcb66680be231e016b Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 19:31:31 +0800 Subject: [PATCH 06/14] update makefile for socket --- src/makefile | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/makefile b/src/makefile index 091d52c7..8a2dca39 100644 --- a/src/makefile +++ b/src/makefile @@ -13,7 +13,10 @@ USE_DP_SUBEIG = 0 # Set USE_FFTW = 1 to use FFTW for fast Fourier transform in vdWDF. Don't open it together with USE_MKL USE_FFTW = 0 # Set DEBUG_MODE = 1 to run with debug mode and print debug output -DEBUG_MODE = 1 +DEBUG_MODE = 0 + +# Set USE_SOCKET = 0 to disable compilation with socket support +USE_SOCKET = 1 # Enable SIMD vectorization for complex stencil routines # CAUTION: for some compilers this results in wrong results! Use for intel/19.0.3 or later versions @@ -109,6 +112,18 @@ OBJSC = main.o initialization.o readfiles.o atomdata.o parallelization.o relax.o mlff/descriptor.o mlff/ddbp_tools.o highT/mlff/internal_energy_model.o cyclix/mlff/cyclix_mlff_tools.o \ mlff/hnl_soap.o +# Switch the socket compilation +# For future developers: please keep the socket options separate from the +# main objects +# Please do not add socket related *.o files manually in OBJSC +ifeq ($(USE_SOCKET), 1) +CPPFLAGS += -Isocket/ -DUSE_SOCKET +OBJSC := $(patsubst main.o,main_socket.o,$(OBJSC)) +OBJSC += socket/driver.o socket/libinetsocket.o socket/libunixsocket.o +endif + + + LIBBASE = ../lib/sparc TESTBASE = ../.ci From 17818956bbe094bab8c28c2e309204f932103b18 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 19:47:44 +0800 Subject: [PATCH 07/14] add debug mode in conda workflow default --- .conda/build.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.conda/build.sh b/.conda/build.sh index 3b689811..23d4b023 100644 --- a/.conda/build.sh +++ b/.conda/build.sh @@ -4,7 +4,7 @@ echo Compiling sparc using ${CPU_COUNT} cores cd ./src make clean -make -j ${CPU_COUNT} USE_MKL=0 USE_SCALAPACK=1 USE_FFTW=1 USE_SOCKET=1 +make -j ${CPU_COUNT} USE_MKL=0 USE_SCALAPACK=1 USE_FFTW=1 USE_SOCKET=1 DEBUG_MODE=1 # ls -al lib echo "Installing sparc into $PREFIX/bin" From dc949798498fcb6940bfb571f8bbea7702c5705b Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 20:12:56 +0800 Subject: [PATCH 08/14] multiple checks for the initialization --- .github/workflows/build-test.yml | 6 +++++- src/makefile | 1 + src/socket/driver.c | 4 +++- 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 81be871f..00cb8623 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -60,5 +60,9 @@ jobs: rm -rf output - name: Test socket compatibility run: | - mpirun -n 1 sparc || true + mpirun -n 1 sparc | tee command-socket.out || true + # Output from the sparc test command + # should contain the -socket subcommand + grep "-socket" command-socket.out + diff --git a/src/makefile b/src/makefile index 8a2dca39..075e9c58 100644 --- a/src/makefile +++ b/src/makefile @@ -141,5 +141,6 @@ sparc: $(OBJSC) .PHONY: clean clean: rm -f $(OBJSC) $(LIBBASE) + rm -f socket/*.o test: ../tests/SPARC_testing_script.py cd ../tests; python SPARC_testing_script.py diff --git a/src/socket/driver.c b/src/socket/driver.c index 95929c1c..2edb074d 100644 --- a/src/socket/driver.c +++ b/src/socket/driver.c @@ -130,8 +130,10 @@ int initialize_Socket(SPARC_OBJ *pSPARC) pSPARC->SocketSCFCount = 0; #ifdef DEBUG - if (rank == 0) + if (rank == 0){ printf("##########################Initializing socket##########################\n"); + printf("SPARC socket_inet type? %d\n", pSPARC->socket_inet); + } #endif // Create a fd according to the socket type int socket_fd = -1; From 890ef39fbdfb52873d8ff21c6198073c3ee1c781 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 20:37:50 +0800 Subject: [PATCH 09/14] add checks for test all --- .github/workflows/build-test.yml | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 00cb8623..513d75ac 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -62,7 +62,13 @@ jobs: run: | mpirun -n 1 sparc | tee command-socket.out || true # Output from the sparc test command - # should contain the -socket subcommand - grep "-socket" command-socket.out + # should contain the -socket subcommand (use regex) + grep "\\-socket" command-socket.out + - name: Run simple socket calculation + run: | + # Requires a local copy of sparc-x-api to function + pip install git+https://github.com/SPARC-X/SPARC-X-API.git@ee3feebfa8a84f72daa52f27e68dbb0d2bcc56db + cd tests/Socket/Al_volume_change + python test_all.py From b6551bda37299b334670b54fe05f156fd7e15a3a Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 20:50:10 +0800 Subject: [PATCH 10/14] use conda to install sparc-x-api --- .github/workflows/build-test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 513d75ac..965fdbce 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -67,7 +67,7 @@ jobs: - name: Run simple socket calculation run: | # Requires a local copy of sparc-x-api to function - pip install git+https://github.com/SPARC-X/SPARC-X-API.git@ee3feebfa8a84f72daa52f27e68dbb0d2bcc56db + mamba install -c conda-forge sparc-x-api==1.0.4 cd tests/Socket/Al_volume_change python test_all.py From 8b0d587f6243d346d6033d53842600a64a619c20 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 20:55:25 +0800 Subject: [PATCH 11/14] remove the quick socket tests --- .github/workflows/build-test.yml | 6 ------ 1 file changed, 6 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 965fdbce..d0fa7031 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -64,11 +64,5 @@ jobs: # Output from the sparc test command # should contain the -socket subcommand (use regex) grep "\\-socket" command-socket.out - - name: Run simple socket calculation - run: | - # Requires a local copy of sparc-x-api to function - mamba install -c conda-forge sparc-x-api==1.0.4 - cd tests/Socket/Al_volume_change - python test_all.py From 9e22c7324de44c2b0d32645a904bb1aa1f4e5f76 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 21:07:16 +0800 Subject: [PATCH 12/14] update changelog --- ChangeLog | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/ChangeLog b/ChangeLog index e212e7b9..ff1c2e58 100644 --- a/ChangeLog +++ b/ChangeLog @@ -2,6 +2,17 @@ -Date -Name -changes +-------------- +Nov 18, 2024 +Name: Tian Tian, Lucas Timmerman +Changes: (initialization.c, initialization.h, electronicGroundState.c, + include/isddft.h, main.c, main_socket.c, makefile, md.c, relax.c, + socket/*.h, socket/*.c, tests/Socket/*, CI workflow, doc) +1. Add socket communication layer and socket submodule +2. Unifying MLFF and DFT single point calculations to Calculate_Properties function +3. Update CI workflow for compiling socket code +4. Update socket mode simple tests +4. Update documentation for socket mode -------------- Nov 13, 2024 From 1d916d14ecaa498e7c49a1abfbf87e6f845d4230 Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Mon, 18 Nov 2024 21:07:54 +0800 Subject: [PATCH 13/14] changelog typo --- ChangeLog | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ChangeLog b/ChangeLog index ff1c2e58..585fab0d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -11,8 +11,8 @@ Changes: (initialization.c, initialization.h, electronicGroundState.c, 1. Add socket communication layer and socket submodule 2. Unifying MLFF and DFT single point calculations to Calculate_Properties function 3. Update CI workflow for compiling socket code -4. Update socket mode simple tests -4. Update documentation for socket mode +4. Update socket mode simple tests +5. Update documentation for socket mode -------------- Nov 13, 2024 From 4c52587ee9bba69da4756ce7abfc6fdb038bd35d Mon Sep 17 00:00:00 2001 From: "T.Tian" Date: Tue, 19 Nov 2024 14:43:46 +0800 Subject: [PATCH 14/14] declaration for calculate_properties in header --- src/include/electronicGroundState.h | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/include/electronicGroundState.h b/src/include/electronicGroundState.h index 1a7a5a3b..a27188f6 100644 --- a/src/include/electronicGroundState.h +++ b/src/include/electronicGroundState.h @@ -15,6 +15,14 @@ #include "isddft.h" +/** + * @brief Calculate properties for fixed atom positions using DFT or MLFF + * + * @param pSPARC + */ +void Calculate_Properties(SPARC_OBJ *pSPARC); + + /** * @brief Calculate the ground state energy and forces for fixed atom positions. *