The sc-config
program can be use to obtain the compilers, compiler options, and libraries needed to use the SC toolkit from your program. This utility is discussed below, along with how the SC toolkit must be initialized in your main
subroutine.
The sc-config program returns information about how SC was compiled and installed. The following information is available:
--prefix
--version
--libdir
--libs
--cppflags
--cc
--cflags
--cxx
--cxxflags
--f77
--f77flags
To use the sc-config program to link your executable to SC, use a Makefile for GNU make similar to the following:
SCCONFIG = /usr/local/mpqc/current/bin/sc-config CXX := $(shell $(SCCONFIG) --cxx) CXXFLAGS := $(shell $(SCCONFIG) --cxxflags) CPPFLAGS := $(shell $(SCCONFIG) --cppflags) LIBS := $(shell $(SCCONFIG) --libs)
myprog: myprog.o $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS)
First the execution environment must be initialized using the ExEnv init member.
ExEnv::init(argc, argv);
By default, all output will go to the console stream, cout. To change this, use the following code:
ostream *outstream = new ofstream(outputfilename); ExEnv::set_out(outstream);
MPI is allowed wait until MPI_Init is called to fill in argc and argv, so you may have to call MPI_Init before you even know that we ready to construct MPIMessageGrp. So if an MPIMessageGrp is needed, it is up to the developer to call MPI_Init to get the argument list for certain MPI implementations.
MPI_Init(&argc, &argv);
When files are read and written, an extension is added to a basename to construct the file name. The default is "SC". To use another basename, make the following call, where basename
is a const char *
:
SCFormIO::set_default_basename(basename);
If your job might run in parallel, then make the following call or the nodes will print redundant information. The myproc
argument is the rank of the called node.
SCFormIO::init_mp(myproc);
This segment of code sets up an object object to provide multi-threading:
RefThreadGrp thread = ThreadGrp::initial_threadgrp(argc, argv); ThreadGrp::set_default_threadgrp(thread); if (thread.nonnull()) ThreadGrp::set_default_threadgrp(thread); else thread = ThreadGrp::get_default_threadgrp();
This segment of code sets up the message passing object:
RefMessageGrp grp = MessageGrp::initial_messagegrp(argc, argv); if (grp.nonnull()) MessageGrp::set_default_messagegrp(grp); else grp = MessageGrp::get_default_messagegrp();
This example code illustrates a complete MP2 energy implementation using the SC Toolkit. First an MP2 class is declared and the necesary base class member functions are provided. Next a ClassDesc is defined. Finally, the member functions are defined.
Note that no main routine is provided. This is because this file is designed to be used to extend the functionality of the mpqc executable. To generate a new mpqc executable with the new class available for use, see the MP2 Implementation Example: Makefile section.
#include <chemistry/qc/wfn/obwfn.h> #include <chemistry/qc/scf/clhf.h> using namespace std; using namespace sc; class MP2: public Wavefunction { Ref<OneBodyWavefunction> ref_mp2_wfn_; double compute_mp2_energy(); public: MP2(const Ref<KeyVal> &); MP2(StateIn &); void save_data_state(StateOut &); void compute(void); void obsolete(void); int nelectron(void); RefSymmSCMatrix density(void); int spin_polarized(void); int value_implemented(void) const; }; static ClassDesc MP2_cd(typeid(MP2), "MP2", 1, "public Wavefunction", 0, create<MP2>, create<MP2>); MP2::MP2(const Ref<KeyVal> &keyval):Wavefunction(keyval) { ref_mp2_wfn_ << keyval->describedclassvalue("reference"); if(ref_mp2_wfn_.null()) { ExEnv::out0() << "reference is null" << endl; abort(); } } MP2::MP2(StateIn &statein):Wavefunction(statein) { ref_mp2_wfn_ << SavableState::restore_state(statein); } void MP2::save_data_state(StateOut &stateout) { Wavefunction::save_data_state(stateout); SavableState::save_state(ref_mp2_wfn_.pointer(),stateout); } void MP2::compute(void) { if(gradient_needed()) { ExEnv::out0() << "No gradients yet" << endl; abort(); } double extra_hf_acc = 10.; ref_mp2_wfn_->set_desired_value_accuracy(desired_value_accuracy() / extra_hf_acc); double refenergy = ref_mp2_wfn_->energy(); double mp2energy = compute_mp2_energy(); ExEnv::out0() << indent << "MP2 Energy = " << mp2energy << endl; set_value(refenergy + mp2energy); set_actual_value_accuracy(ref_mp2_wfn_->actual_value_accuracy() * extra_hf_acc); } void MP2::obsolete(void) { Wavefunction::obsolete(); ref_mp2_wfn_->obsolete(); } int MP2::nelectron(void) { return ref_mp2_wfn_->nelectron(); } RefSymmSCMatrix MP2::density(void) { ExEnv::out0() << "No density yet" << endl; abort(); return 0; } int MP2::spin_polarized(void) { return 0; } int MP2::value_implemented(void) const { return 1; } double MP2::compute_mp2_energy() { if(molecule()->point_group()->char_table().order() != 1) { ExEnv::out0() << "C1 symmetry only" << endl; abort(); } RefSCMatrix vec = ref_mp2_wfn_->eigenvectors(); int nao = vec.nrow(); int nmo = vec.ncol(); int nocc = ref_mp2_wfn_->nelectron()/2; int nvir = nmo - nocc; double *cvec = new double [vec.nrow() * vec.ncol()]; vec->convert(cvec); double *pqrs = new double [nao * nao * nao * nao]; for(int n = 0; n < nao*nao*nao*nao; n++) pqrs[n] = 0.0; Ref<TwoBodyInt> twoint = integral()->electron_repulsion(); const double *buffer = twoint->buffer(); Ref<GaussianBasisSet> basis = this->basis(); int nshell = basis->nshell(); for(int P = 0; P < nshell; P++) { int nump = basis->shell(P).nfunction(); for(int Q = 0; Q < nshell; Q++) { int numq = basis->shell(Q).nfunction(); for(int R = 0; R < nshell; R++) { int numr = basis->shell(R).nfunction(); for(int S = 0; S < nshell; S++) { int nums = basis->shell(S).nfunction(); twoint->compute_shell(P,Q,R,S); int index = 0; for(int p=0; p < nump; p++) { int op = basis->shell_to_function(P)+p; for(int q = 0; q < numq; q++) { int oq = basis->shell_to_function(Q)+q; for(int r = 0; r < numr; r++) { int oor = basis->shell_to_function(R)+r; for(int s = 0; s < nums; s++,index++) { int os = basis->shell_to_function(S)+s; int ipqrs = (((op*nao+oq)*nao+oor)*nao+os); pqrs[ipqrs] = buffer[index]; } } } } } } } } twoint = 0; double *ijkl = new double [nmo * nmo * nmo * nmo]; int idx = 0; for(int i = 0; i < nmo; i++) { for(int j = 0; j < nmo; j++) { for(int k = 0; k < nmo; k++) { for(int l = 0; l < nmo; l++, idx++) { ijkl[idx] = 0.0; int index = 0; for(int p = 0; p < nao; p++) { for(int q = 0; q < nao; q++) { for(int r = 0; r < nao; r++) { for(int s = 0; s < nao; s++,index++) { ijkl[idx] += cvec[p*nmo + i] * cvec[q*nmo +j] * cvec[r*nmo + k] * cvec[s*nmo + l] * pqrs[index]; } } } } } } } } delete [] pqrs; delete [] cvec; double *evals = new double [nmo]; ref_mp2_wfn_->eigenvalues()->convert(evals); double energy = 0.0; for(int i=0; i < nocc; i++) { for(int j=0; j < nocc; j++) { for(int a=nocc; a < nmo; a++) { for(int b=nocc; b < nmo; b++) { int iajb = (((i*nmo+a)*nmo+j)*nmo+b); int ibja = (((i*nmo+b)*nmo+j)*nmo+a); energy += (2 * ijkl[iajb] - ijkl[ibja]) * ijkl[iajb]/ (evals[i] + evals[j] - evals[a] - evals[b]); } } } } delete [] ijkl; delete [] evals; return energy; }
00001 00002 #include <chemistry/qc/wfn/obwfn.h> 00003 #include <chemistry/qc/scf/clhf.h> 00004 00005 using namespace std; 00006 using namespace sc; 00007 00008 class MP2: public Wavefunction { 00009 Ref<OneBodyWavefunction> ref_mp2_wfn_; 00010 double compute_mp2_energy(); 00011 00012 public: 00013 MP2(const Ref<KeyVal> &); 00014 MP2(StateIn &); 00015 void save_data_state(StateOut &); 00016 void compute(void); 00017 void obsolete(void); 00018 int nelectron(void); 00019 RefSymmSCMatrix density(void); 00020 int spin_polarized(void); 00021 int value_implemented(void) const; 00022 }; 00023 00024 static ClassDesc MP2_cd(typeid(MP2), "MP2", 1, "public Wavefunction", 00025 0, create<MP2>, create<MP2>); 00026 00027 MP2::MP2(const Ref<KeyVal> &keyval):Wavefunction(keyval) { 00028 ref_mp2_wfn_ << keyval->describedclassvalue("reference"); 00029 if(ref_mp2_wfn_.null()) { 00030 ExEnv::out0() << "reference is null" << endl; 00031 abort(); 00032 } 00033 } 00034 00035 MP2::MP2(StateIn &statein):Wavefunction(statein) 00036 { 00037 ref_mp2_wfn_ << SavableState::restore_state(statein); 00038 } 00039 00040 void 00041 MP2::save_data_state(StateOut &stateout) { 00042 Wavefunction::save_data_state(stateout); 00043 00044 SavableState::save_state(ref_mp2_wfn_.pointer(),stateout); 00045 } 00046 00047 void 00048 MP2::compute(void) 00049 { 00050 if(gradient_needed()) { 00051 ExEnv::out0() << "No gradients yet" << endl; 00052 abort(); 00053 } 00054 00055 double extra_hf_acc = 10.; 00056 ref_mp2_wfn_->set_desired_value_accuracy(desired_value_accuracy() 00057 / extra_hf_acc); 00058 double refenergy = ref_mp2_wfn_->energy(); 00059 double mp2energy = compute_mp2_energy(); 00060 00061 ExEnv::out0() << indent << "MP2 Energy = " << mp2energy << endl; 00062 00063 set_value(refenergy + mp2energy); 00064 set_actual_value_accuracy(ref_mp2_wfn_->actual_value_accuracy() 00065 * extra_hf_acc); 00066 } 00067 00068 void 00069 MP2::obsolete(void) { 00070 Wavefunction::obsolete(); 00071 ref_mp2_wfn_->obsolete(); 00072 } 00073 00074 int 00075 MP2::nelectron(void) { 00076 return ref_mp2_wfn_->nelectron(); 00077 } 00078 00079 RefSymmSCMatrix 00080 MP2::density(void) { 00081 ExEnv::out0() << "No density yet" << endl; 00082 abort(); 00083 return 0; 00084 } 00085 00086 int 00087 MP2::spin_polarized(void) { 00088 return 0; 00089 } 00090 00091 int 00092 MP2::value_implemented(void) const { 00093 return 1; 00094 } 00095 00096 double 00097 MP2::compute_mp2_energy() 00098 { 00099 if(molecule()->point_group()->char_table().order() != 1) { 00100 ExEnv::out0() << "C1 symmetry only" << endl; 00101 abort(); 00102 } 00103 00104 RefSCMatrix vec = ref_mp2_wfn_->eigenvectors(); 00105 00106 int nao = vec.nrow(); 00107 int nmo = vec.ncol(); 00108 int nocc = ref_mp2_wfn_->nelectron()/2; 00109 int nvir = nmo - nocc; 00110 00111 double *cvec = new double [vec.nrow() * vec.ncol()]; 00112 vec->convert(cvec); 00113 00114 double *pqrs = new double [nao * nao * nao * nao]; 00115 for(int n = 0; n < nao*nao*nao*nao; n++) pqrs[n] = 0.0; 00116 00117 Ref<TwoBodyInt> twoint = integral()->electron_repulsion(); 00118 const double *buffer = twoint->buffer(); 00119 00120 Ref<GaussianBasisSet> basis = this->basis(); 00121 00122 int nshell = basis->nshell(); 00123 for(int P = 0; P < nshell; P++) { 00124 int nump = basis->shell(P).nfunction(); 00125 00126 for(int Q = 0; Q < nshell; Q++) { 00127 int numq = basis->shell(Q).nfunction(); 00128 00129 for(int R = 0; R < nshell; R++) { 00130 int numr = basis->shell(R).nfunction(); 00131 00132 for(int S = 0; S < nshell; S++) { 00133 int nums = basis->shell(S).nfunction(); 00134 00135 twoint->compute_shell(P,Q,R,S); 00136 00137 int index = 0; 00138 for(int p=0; p < nump; p++) { 00139 int op = basis->shell_to_function(P)+p; 00140 00141 for(int q = 0; q < numq; q++) { 00142 int oq = basis->shell_to_function(Q)+q; 00143 00144 for(int r = 0; r < numr; r++) { 00145 int oor = basis->shell_to_function(R)+r; 00146 00147 for(int s = 0; s < nums; s++,index++) { 00148 int os = basis->shell_to_function(S)+s; 00149 00150 int ipqrs = (((op*nao+oq)*nao+oor)*nao+os); 00151 00152 pqrs[ipqrs] = buffer[index]; 00153 00154 } 00155 } 00156 } 00157 } 00158 00159 } 00160 } 00161 } 00162 } 00163 00164 twoint = 0; 00165 00166 double *ijkl = new double [nmo * nmo * nmo * nmo]; 00167 00168 int idx = 0; 00169 for(int i = 0; i < nmo; i++) { 00170 for(int j = 0; j < nmo; j++) { 00171 for(int k = 0; k < nmo; k++) { 00172 for(int l = 0; l < nmo; l++, idx++) { 00173 00174 ijkl[idx] = 0.0; 00175 00176 int index = 0; 00177 for(int p = 0; p < nao; p++) { 00178 for(int q = 0; q < nao; q++) { 00179 for(int r = 0; r < nao; r++) { 00180 for(int s = 0; s < nao; s++,index++) { 00181 00182 ijkl[idx] += cvec[p*nmo + i] * cvec[q*nmo +j] 00183 * cvec[r*nmo + k] * cvec[s*nmo + l] 00184 * pqrs[index]; 00185 00186 } 00187 } 00188 } 00189 } 00190 00191 } 00192 } 00193 } 00194 } 00195 00196 delete [] pqrs; 00197 delete [] cvec; 00198 00199 double *evals = new double [nmo]; 00200 ref_mp2_wfn_->eigenvalues()->convert(evals); 00201 00202 double energy = 0.0; 00203 for(int i=0; i < nocc; i++) { 00204 for(int j=0; j < nocc; j++) { 00205 for(int a=nocc; a < nmo; a++) { 00206 for(int b=nocc; b < nmo; b++) { 00207 00208 int iajb = (((i*nmo+a)*nmo+j)*nmo+b); 00209 int ibja = (((i*nmo+b)*nmo+j)*nmo+a); 00210 00211 energy += (2 * ijkl[iajb] - ijkl[ibja]) * ijkl[iajb]/ 00212 (evals[i] + evals[j] - evals[a] - evals[b]); 00213 00214 } 00215 } 00216 } 00217 } 00218 00219 delete [] ijkl; 00220 delete [] evals; 00221 00222 return energy; 00223 }
This example Makefile demonstrates how to link in a new class to form a new mpqc executable, here named mp2. The code is given in the MP2 Implementation Example: Source section. The sc-config command is used to obtain information about how the SC toolkit was compiled and installed. The library specified with -lmpqc provides the main routine from mpqc.
# Change this to the path to your installed sc-config script. SCCONFIG = /usr/local/mpqc/current/bin/sc-config CXX := $(shell $(SCCONFIG) --cxx) CXXFLAGS := $(shell $(SCCONFIG) --cxxflags) CPPFLAGS := $(shell $(SCCONFIG) --cppflags) LIBS := $(shell $(SCCONFIG) --libs) LIBDIR := $(shell $(SCCONFIG) --libdir) mp2: mp2.o $(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) -lmpqc $(LIBS)
00001 # Change this to the path to your installed sc-config script. 00002 SCCONFIG = /usr/local/mpqc/current/bin/sc-config 00003 CXX := $(shell $(SCCONFIG) --cxx) 00004 CXXFLAGS := $(shell $(SCCONFIG) --cxxflags) 00005 CPPFLAGS := $(shell $(SCCONFIG) --cppflags) 00006 LIBS := $(shell $(SCCONFIG) --libs) 00007 LIBDIR := $(shell $(SCCONFIG) --libdir) 00008 00009 mp2: mp2.o 00010 $(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) -lmpqc $(LIBS)
This input file can be used with the program illustrated in the MP2 Implementation Example: Source section. It will compute the MP2 energy using the new MP2 class. Note that only the object-oriented input format can be used with user provided classes.
% emacs should use -*- KeyVal -*- mode molecule<Molecule>: ( symmetry = C1 unit = angstrom { atoms geometry } = { O [ 0.00000000 0.00000000 0.37000000 ] H [ 0.78000000 0.00000000 -0.18000000 ] H [ -0.78000000 0.00000000 -0.18000000 ] } ) basis<GaussianBasisSet>: ( name = "STO-3G" molecule = $:molecule ) mpqc: ( checkpoint = no savestate = no % MP2 is the new class. Change MP2 to MBPT2 to test % against the standard MP2 code mole<MP2>: ( molecule = $:molecule basis = $:basis reference<CLHF>: ( molecule = $:molecule basis = $:basis memory = 16000000 ) ) )
00001 % emacs should use -*- KeyVal -*- mode 00002 molecule<Molecule>: ( 00003 symmetry = C1 00004 unit = angstrom 00005 { atoms geometry } = { 00006 O [ 0.00000000 0.00000000 0.37000000 ] 00007 H [ 0.78000000 0.00000000 -0.18000000 ] 00008 H [ -0.78000000 0.00000000 -0.18000000 ] 00009 } 00010 ) 00011 basis<GaussianBasisSet>: ( 00012 name = "STO-3G" 00013 molecule = $:molecule 00014 ) 00015 mpqc: ( 00016 checkpoint = no 00017 savestate = no 00018 % MP2 is the new class. Change MP2 to MBPT2 to test 00019 % against the standard MP2 code 00020 mole<MP2>: ( 00021 molecule = $:molecule 00022 basis = $:basis 00023 reference<CLHF>: ( 00024 molecule = $:molecule 00025 basis = $:basis 00026 memory = 16000000 00027 ) 00028 ) 00029 ) 00030