Commit eec45811 authored by Xiaolin Liu's avatar Xiaolin Liu Committed by Andrew Wildman
Browse files

Implement CCSD

parent 12c08895
......@@ -16,6 +16,8 @@ external/libxc/lib
external/libxc/share
external/libxc/src
external/libxc/tmp
external/tiledarray/
doc
*.*.swp
*~
......@@ -12,7 +12,7 @@ build_intel:
- source /home/ci_software/intel-19.1.sh
- mkdir build
- cd build
- cmake -DLibint2_ROOT=/sw/libint/2/5/0/intel/19/1 =-DCMAKE_CXX_FLAGS='-O3' -DCMAKE_C_FLAGS='-O3' -DCMAKE_Fortran_FLAGS='-O3' ..
- cmake -DLibint2_ROOT=/sw/libint/2/5/0/intel/19/1 =-DCMAKE_CXX_FLAGS='-O3' -DCMAKE_C_FLAGS='-O3' -DCMAKE_Fortran_FLAGS='-O3' ..
- make -j5
- ctest
......@@ -42,6 +42,6 @@ build_mpi_gcc:
- source /home/ci_software/gcc-8.2.0_mpich-3.2.1.sh
- mkdir build
- cd build
- cmake -DLibint2_ROOT=/sw/libint/2/5/0/gcc/8/2/0 -DCMAKE_CXX_FLAGS='-O3' -DCMAKE_C_FLAGS='-O3' -DCMAKE_Fortran_FLAGS='-O3' -DCQ_ENABLE_MPI=ON -DCQ_SCALAPACK_LIBRARIES="$SCALAPACK_LIBRARY_DIR/libscalapack.a" -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release ..
- cmake -DLibint2_ROOT=/sw/libint/2/5/0/gcc/8/2/0 -DCMAKE_CXX_FLAGS='-O3' -DCMAKE_C_FLAGS='-O3' -DCMAKE_Fortran_FLAGS='-O3' -DCQ_ENABLE_MPI=ON -DCQ_SCALAPACK_LIBRARIES="$SCALAPACK_LIBRARY_DIR/libscalapack.a" -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release ..
- make -j6
- make test
......@@ -41,6 +41,7 @@ option(CQ_ENABLE_MPI "Enable MPI parallelism" OFF)
option(ENABLE_COVERAGE "Enable coverage and profiling" OFF)
option(CQ_LINALG_USESYSTEM "Use System defaults for LA libs" OFF)
option(CQ_EXTERNAL_OPENMP "Force linking to an external OpenMP library" OFF)
option(CQ_ENABLE_TA "Enable TA" OFF)
# Libint building options
......@@ -159,7 +160,16 @@ list(APPEND CQEX_LINK ChronusQ::HDF5 )
message( "\n\n" )
if( CQ_ENABLE_TA )
if(NOT CQ_ENABLE_MPI)
message(FATAL_ERROR "TiledArray must be compiled with mpi!")
endif()
message("Enable TiledArray")
include(FindTiledArray)
set( CQ_HAS_TA ON )
else()
message("TiledArray will not be compiled, Coupled-Cluster will be disabled.")
endif()
# Compiling with coverage report
if(ENABLE_COVERAGE)
......
#
# This file is part of the Chronus Quantum (ChronusQ) software package
#
# Copyright (C) 2014-2020 Li Research Group (University of Washington)
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Contact the Developers:
# E-Mail: xsli@uw.edu
#
#
#only support openblas built by CQ now.
if(NOT CQ_NEED_OPENBLAS)
message ( FATAL_ERROR "Only OPENBLAS is supported!" )
endif()
include(ExternalProject)
set( TA_PREFIX ${PROJECT_SOURCE_DIR}/external/tiledarray )
set( TA_INCLUDEDIR ${TA_PREFIX}/include )
set( TA_LIBDIR ${TA_PREFIX}/lib )
include_directories("${TA_INCLUDEDIR}")
link_directories("${TA_LIBDIR}")
if( NOT EXISTS "${TA_INCLUDEDIR}/tiledarray.h" )
ExternalProject_Add(tiledarray
PREFIX ${TA_PREFIX}
GIT_REPOSITORY https://github.com/ValeevGroup/tiledarray.git
GIT_TAG 2aef4e6552a9ae5f4b1a4dbdaf1346adf3610abd
CMAKE_ARGS
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-DMPI_C_COMPILER=${MPI_C_COMPILER}
-DMPI_CXX_COMPILER=${MPI_CXX_COMPILER}
-DCMAKE_INSTALL_PREFIX=${TA_PREFIX}
-DEIGEN3_INCLUDE_DIR=${EIGEN3_INCLUDE_DIR}
-DENABLE_MKL=OFF
-D LAPACK_LIBRARIES="${PROJECT_SOURCE_DIR}/external/openblas/lib/libopenblas.a -lm -lgfortran"
BUILD_IN_SOURCE 1
UPDATE_COMMAND echo 'NO UPDATE'
)
if( TARGET openblas )
add_dependencies( tiledarray openblas )
endif()
list(APPEND CQEX_DEP tiledarray)
message( STATUS "Opting to build a copy of TiledArray")
else()
message( STATUS "Found TiledArray")
endif()
list(APPEND CQEX_LINK ${TA_LIBDIR}/libtiledarray.so)
list(APPEND CQEX_LINK ${TA_LIBDIR}/libMADworld.so)
......@@ -46,7 +46,7 @@
// MPI / ScaLAPACK / BLACS
#cmakedefine CQ_ENABLE_MPI
#cmakedefine CQ_HAS_TA
/*
// Add a fatal preprocessor error if less than ICC 17.0.4
#ifdef __INTEL_COMPILER
......
/*
* This file is part of the Chronus Quantum (ChronusQ) software package
*
* Copyright (C) 2014-2020 Li Research Group (University of Washington)
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Contact the Developers:
* E-Mail: xsli@uw.edu
*
*/
#ifdef CQ_HAS_TA
#pragma once
#include <singleslater.hpp>
#include <electronintegrals.hpp>
#include <wavefunction.hpp>
#include <chronusq_sys.hpp>
#include <singleslater/hartreefock.hpp>
#include <util/files.hpp>
#include <util/math.hpp>
#include <integrals.hpp>
#include <cerr.hpp>
#include <tiledarray.h>
#include "./coupledcluster/DiskDIIS.hpp"
namespace ChronusQ {
enum class CC_TYPE { CCSD};
struct CoupledClusterSettings {
CC_TYPE cctype = CC_TYPE::CCSD;
double eConv = 1e-8;
double tConv = 1e-6;
int maxiter = 1000;
};
struct CCBase
{
virtual void getCorrEnergy() = 0;
virtual void run() = 0;
bool useDIIS = true;
size_t nDIIS = 8;
CoupledClusterSettings ccSettings;
};
template <typename MatsT, typename IntsT>
class CCSD : public CCBase
{
using TArray = TA::TArray<MatsT>;
protected:
TArray T1_;
TArray T2_;
SingleSlater<MatsT,IntsT>& ref_;
public:
CCSD(SingleSlater<MatsT,IntsT>& ref): ref_(ref){}
~CCSD(){};
std::map<std::string,TArray> antiSymMoints;
std::map<std::string,TArray> fockMatrix_ta;
double CorrE = 0.0;
void getCorrEnergy() override;
void updateT1(const TArray T1_old, const TArray T2_old);
void updateT2(const TArray T1_old, TArray T2_old);
void transformInts();
void run() override ;
void runConventional();
//Auxiliary variables and functions to help initialization of TA tensors
TA::TiledRange1 orange;
TA::TiledRange1 vrange;
std::size_t blksize{4};
void initRanges();
//Singles intermediates and tensors
void formA_1();
void formA_2();
void formA_3();
void formA_4();
TArray A_1;
TArray A_2;
TArray A_3;
TArray A_4;
//Doubles intermediates and tensors
void formB_1();
void formB_2();
void formB_3();
void formB_5();
void formB_6();
TArray B_1;
TArray B_2;
TArray B_3;
TArray B_5;
TArray B_6;
void initIntermediates();
void initAmplitudes();
void initFock();
void doDIIS(size_t t1offset, size_t NT1blks, size_t NT2blks, TArray T1_old, TArray T2_old, std::shared_ptr<DiskDIIS<MatsT>> diis );
};// namespace CCSD
}; // namespace ChronusQ
#endif
This diff is collapsed.
/*
* This file is part of the Chronus Quantum (ChronusQ) software package
*
* Copyright (C) 2014-2020 Li Research Group (University of Washington)
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Contact the Developers:
* E-Mail: xsli@uw.edu
*
*/
#pragma once
#include <cqlinalg/blas1.hpp>
#include <cqlinalg/solve.hpp>
#include <util/files.hpp>
#include <memmanager.hpp>
namespace ChronusQ {
/**
* \brief A general DIIS class for extrapolations
* based on a series of error metrics stored on disk.
*
*/
template <typename T>
class DiskDIIS{
public:
/**
* DiskDIIS Constructor. Constructs a DiskDIIS object
*
* \param [in] n Dimension of extrapolation space. Two temporary arrays of length n will be created
* \param [in] maxdiis Maximum number of error vectors
* \param [in] savFile File manager
* \param [in] memManager Memory manager
*/
DiskDIIS(int n, int maxdiis, SafeFile savFile, CQMemManager & memManager);
~DiskDIIS();
/**
* getVector. Get current solution
*
* \param [out] vector Pointer to target array
* \param [in] dim Number of elements to be copied to vector
* \param [in] off Offset in full array (if multiple quantities extrapolated simultaneously)
*/
void getVector(T * vector,int dim, int off);
/**
* setVector. Set current solution
*
* \param [in] vector Pointer to current solution
* \param [in] dim Number of elements to be copied to DiskDIIS
* \param [in] off Offset in full array (if multiple quantities extrapolated simultaneously)
*/
void setVector(T * vector,int dim, int off);
/**
* setErrorVector. Set current error vector
*
* \param [in] vector Pointer to current error vector
* \param [in] dim Number of elements to be copied to DiskDIIS
* \param [in] off Offset in full array (if multiple quantities extrapolated simultaneously)
*/
void setErrorVector(T * vector,int dim, int off);
/**
* extrapolate. solve DiskDIIS equations and extrapolation solution
*
* \param [in] vector Pointer to current error vector
* \param [in] dim Number of elements to be copied to DiskDIIS
* \param [in] off Offset in full array (if multiple quantities extrapolated simultaneously)
*/
void extrapolate();
protected:
/// file manager
SafeFile savFile_;
/// memory manager
CQMemManager & memManager_;
/// write current error vector to disk
void write_error_vector(T * vector);
/// write current solution to disk
void write_vector(T * vector);
/// DiskDIIS temporary storage 1
T * tmp1_;
/// DiskDIIS temporary storage 2
T * tmp2_;
/// determine DiskDIIS expansion coefficients
void DIISCoefficients(int nvec);
/// maximum number of vectors
int maxdiis_;
/// array containing DiskDIIS expansion coefficients
T * diisvec_;
/// dimension of the solution/error vector
int dimdiis_;
/// current number of DiskDIIS vectors
int diis_iter_;
/// DiskDIIS vector to be replaced
int replace_diis_iter_;
};
template <typename T>
DiskDIIS<T>::DiskDIIS(int n, int maxdiis, SafeFile savFile, CQMemManager & memManager):
savFile_(savFile), memManager_(memManager)
{
dimdiis_ = n;
maxdiis_ = maxdiis;
diis_iter_ = 0;
replace_diis_iter_ = 1;
diisvec_ = memManager_.template malloc<T>(maxdiis_+1);
tmp1_ = memManager_.template malloc<T>(dimdiis_);
tmp2_ = memManager_.template malloc<T>(dimdiis_);
}
template <typename T>
DiskDIIS<T>::~DiskDIIS(){
memManager_.free(diisvec_);
memManager_.free(tmp1_);
memManager_.free(tmp2_);
}
// Write current solution vector to disk
template <typename T>
void DiskDIIS<T>::write_vector(T * vector){
// Name the entry in according to the current DiskDIIS iteration.
// If we already have maxdiis_ vectors, then replace one.
char * oldvector = memManager_.template malloc<char>(1000);
if ( diis_iter_ <= maxdiis_ ){
sprintf(oldvector,"/DIIS/VECTOR%i",diis_iter_);
}
else{
sprintf(oldvector,"/DIIS/VECTOR%i",replace_diis_iter_);
}
// Write the current solution vector.
savFile_.safeWriteData(oldvector, vector, {dimdiis_});
memManager_.free(oldvector);
}
template <typename T>
void DiskDIIS<T>::getVector(T * vector,int dim, int off) {
std::copy_n(tmp1_+off,dim,vector);
}
template <typename T>
void DiskDIIS<T>::setVector(T * vector,int dim, int off) {
std::copy_n(vector,dim,tmp1_+off);
}
template <typename T>
void DiskDIIS<T>::setErrorVector(T * vector,int dim, int off) {
std::copy_n(vector,dim,tmp2_+off);
}
// Write current error vector to disk
template <typename T>
void DiskDIIS<T>::write_error_vector(T * vector){
// Name the entry in according to the current DiskDIIS iteration.
// If we already have maxdiis_ vectors, then replace one.
char * evector = memManager_.template malloc<char>(1000);
if ( diis_iter_ <= maxdiis_ ){
sprintf(evector,"/DIIS/ERROR%i",diis_iter_);
}
else{
sprintf(evector,"/DIIS/ERROR%i",replace_diis_iter_);
}
if ( diis_iter_ == 0 ) {
// On the first iteration, write an entry that will hold the error matrix
T * temp = (T*)malloc(maxdiis_*maxdiis_*sizeof(T));
memset((void*)temp,'\0',maxdiis_*maxdiis_*sizeof(T));
savFile_.safeWriteData("/DIIS/ERROR_MATRIX", temp, {maxdiis_*maxdiis_});
free(temp);
}
// Write current error vector.
savFile_.safeWriteData(evector, vector, {dimdiis_});
memManager_.free(evector);
}
// Perform DIIS extrapolation. solution in tmp1_
template <typename T>
void DiskDIIS<T>::extrapolate() {
write_vector(tmp1_);
write_error_vector(tmp2_);
if ( diis_iter_ > 1 ) {
// Compute coefficients for the extrapolation
DIISCoefficients( diis_iter_ < maxdiis_ ? diis_iter_ : maxdiis_ );
memset((void*)tmp1_,'\0',dimdiis_*sizeof(T));
char * oldvector = memManager_.template malloc<char>(1000);
int max = diis_iter_;
if (max > maxdiis_) max = maxdiis_;
// Read each of the old vectors from disk.
for (int j = 1; j <= max; j++){
sprintf(oldvector,"/DIIS/VECTOR%i",j);
savFile_.readData(oldvector,tmp2_);
// Accumulate extrapolated vector.
AXPY<T>(dimdiis_,diisvec_[j-1],tmp2_,1,tmp1_,1);
}
memManager_.free(oldvector);
}
if (diis_iter_ <= maxdiis_){
diis_iter_++;
}
//else {
// // If we already have maxdiis_ vectors, choose the one with
// // the largest error as the one to replace.
// int jmax = 1;
// double max = -1.0e99;
// char * evector = memManager_.template malloc<char>(1000);
// for (int j = 1; j <= maxdiis_; j++){
// sprintf(evector,"/DIIS/ERROR%i",j);
// savFile_.readData(evector,tmp2_);
// double nrm = sqrt(std::real(InnerProd<T>(dimdiis_,tmp2_,1,tmp2_,1))); // couldn't get TwoNorm template to work...
// if ( nrm > max ) {
// max = nrm;
// jmax = j;
// }
// }
// replace_diis_iter_ = jmax;
// memManager_.free(evector);
//}
else if (replace_diis_iter_ < maxdiis_) {
replace_diis_iter_++;
}else {
replace_diis_iter_ = 1;
}
}
// Evaluate extrapolation coefficients for DiskDIIS.
template <typename T>
void DiskDIIS<T>::DIISCoefficients(int nvec){
// Allocate memory for small matrices/vectors.
int * ipiv = memManager_.template malloc<int>(nvec+1);
T * temp = memManager_.template malloc<T>(maxdiis_*maxdiis_);
T * A = memManager_.template malloc<T>((nvec+1)*(nvec+1));
T * B = memManager_.template malloc<T>(nvec+1);
memset((void*)A,'\0',(nvec+1)*(nvec+1)*sizeof(T));
memset((void*)B,'\0',(nvec+1)*sizeof(T));
B[nvec] = -1.0;
char * evector = memManager_.template malloc<char>(1000);
// Read in the previous error matrix, so we don't have
// to build the entire thing each iteration.
savFile_.readData("/DIIS/ERROR_MATRIX",temp);
// Reshape the error matrix, in case its dimension is less than maxdiis_.
for (int i = 0; i < nvec; i++){
for (int j = 0; j < nvec; j++){
A[i*(nvec+1)+j] = temp[i*maxdiis_+j];
}
}
if (nvec <= 3) {
// At early iterations, just build the whole matrix.
for (int i = 0; i < nvec; i++) {
sprintf(evector,"/DIIS/ERROR%i",i+1);
savFile_.readData(evector,tmp1_);
for (int j = i+1; j < nvec; j++){
sprintf(evector,"/DIIS/ERROR%i",j+1);
savFile_.readData(evector,tmp2_);
T sum = InnerProd<T>(dimdiis_,tmp1_,1,tmp2_,1);
A[i*(nvec+1)+j] = sum;
A[j*(nvec+1)+i] = sum;
}
T sum = InnerProd<T>(dimdiis_,tmp1_,1,tmp1_,1);
A[i*(nvec+1)+i] = sum;
}
}else {
// At later iterations, don't build the whote matrix.
// Just replace one row/column.
// Which row/column will be replaced?
int i = nvec < maxdiis_ ? nvec - 1 : replace_diis_iter_ - 1;
sprintf(evector,"/DIIS/ERROR%i",i+1);
savFile_.readData(evector,tmp1_);
for (int j = 0; j < nvec; j++){
sprintf(evector,"/DIIS/ERROR%i",j+1);
savFile_.readData(evector,tmp2_);
T sum = InnerProd<T>(dimdiis_,tmp1_,1,tmp2_,1);
A[i*(nvec+1)+j] = sum;
A[j*(nvec+1)+i] = sum;
}
}
int j = nvec;
for (int i = 0; i < (nvec+1); i++){
A[j*(nvec+1)+i] = -1.0;
A[i*(nvec+1)+j] = -1.0;
}
A[(nvec+1)*