~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric : revision 4

src/lib/chemistry/qc/mbptr12/r12int_eval.cc

// r12int_eval.cc

// Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>

// Maintainer: EV

// This file is part of the SC Toolkit.

// The SC Toolkit is free software; you can redistribute it and/or modify

// it under the terms of the GNU Library General Public License as published by

// the Free Software Foundation; either version 2, or (at your option)

// any later version.

// The SC Toolkit 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 Library General Public License for more details.

// You should have received a copy of the GNU Library General Public License

// along with the SC Toolkit; see the file COPYING.LIB. If not, write to

// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.

// The U.S. Government is granted a limited license as per AL 91-7.

#ifdef __GNUG__

#pragma implementation

#endif

#include <util/misc/formio.h>

#include <util/ref/ref.h>

#include <util/state/state_bin.h>

#include <math/scmat/local.h>

#include <chemistry/qc/mbptr12/vxb_eval_info.h>

#include <chemistry/qc/mbptr12/pairiter.h>

#include <chemistry/qc/mbptr12/r12int_eval.h>

using namespace std;

using namespace sc;

#define TEST_FOCK 0

#define NOT_INCLUDE_DIAGONAL_VXB_CONTIBUTIONS 0

inline int max(int a,int b) { return (a > b) ? a : b;}

/*-----------------

R12IntEval

-----------------*/

static ClassDesc R12IntEval_cd(

typeid(R12IntEval),"R12IntEval",1,"virtual public SavableState",

0, 0, 0);

R12IntEval::R12IntEval(const Ref<R12IntEvalInfo>& r12info, bool gbc, bool ebc,

LinearR12::ABSMethod abs_method,

LinearR12::StandardApproximation stdapprox, bool follow_ks_ebcfree) :

r12info_(r12info), gbc_(gbc), ebc_(ebc), abs_method_(abs_method),

stdapprox_(stdapprox), spinadapted_(false), include_mp1_(false),

follow_ks_ebcfree_(follow_ks_ebcfree), debug_(0), evaluated_(false)

{

int nocc_act = r12info_->nocc_act();

int nvir_act = r12info_->nvir_act();

dim_ij_aa_ = new SCDimension((nocc_act*(nocc_act-1))/2);

dim_ij_ab_ = new SCDimension(nocc_act*nocc_act);

dim_ij_s_ = new SCDimension((nocc_act*(nocc_act+1))/2);

dim_ij_t_ = dim_ij_aa_;

dim_ab_aa_ = new SCDimension((nvir_act*(nvir_act-1))/2);

dim_ab_ab_ = new SCDimension(nvir_act*nvir_act);

Ref<LocalSCMatrixKit> local_matrix_kit = new LocalSCMatrixKit();

Vaa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ij_aa_);

Vab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ij_ab_);

Xaa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ij_aa_);

Xab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ij_ab_);

Baa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ij_aa_);

Bab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ij_ab_);

if (ebc_ == false) {

Aaa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

Aab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

if (follow_ks_ebcfree_) {

Ac_aa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

Ac_ab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

}

T2aa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

T2ab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

Raa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

Rab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

}

emp2pair_aa_ = local_matrix_kit->vector(dim_ij_aa_);

emp2pair_ab_ = local_matrix_kit->vector(dim_ij_ab_);

init_intermeds_();

init_tforms_();

}

R12IntEval::R12IntEval(StateIn& si) : SavableState(si)

{

int gbc; si.get(gbc); gbc_ = (bool) gbc;

100

int ebc; si.get(ebc); ebc_ = (bool) ebc;

101

int absmethod; si.get(absmethod); abs_method_ = (LinearR12::ABSMethod) absmethod;

102

int stdapprox; si.get(stdapprox); stdapprox_ = (LinearR12::StandardApproximation) stdapprox;

103

int follow_ks_ebcfree; si.get(follow_ks_ebcfree); follow_ks_ebcfree_ = static_cast<bool>(follow_ks_ebcfree);

104

105

r12info_ << SavableState::restore_state(si);

106

dim_ij_aa_ << SavableState::restore_state(si);

107

dim_ij_ab_ << SavableState::restore_state(si);

108

dim_ij_s_ << SavableState::restore_state(si);

109

dim_ij_t_ << SavableState::restore_state(si);

110

dim_ab_aa_ << SavableState::restore_state(si);

111

dim_ab_ab_ << SavableState::restore_state(si);

112

113

Ref<LocalSCMatrixKit> local_matrix_kit = new LocalSCMatrixKit();

114

Vaa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ij_aa_);

115

Vab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ij_ab_);

116

Xaa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ij_aa_);

117

Xab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ij_ab_);

118

Baa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ij_aa_);

119

Bab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ij_ab_);

120

if (ebc_ == false) {

121

Aaa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

122

Aab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

123

if (follow_ks_ebcfree_) {

124

Ac_aa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

125

Ac_ab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

126

}

127

T2aa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

128

T2ab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

129

Raa_ = local_matrix_kit->matrix(dim_ij_aa_,dim_ab_aa_);

130

Rab_ = local_matrix_kit->matrix(dim_ij_ab_,dim_ab_ab_);

131

}

132

emp2pair_aa_ = local_matrix_kit->vector(dim_ij_aa_);

133

emp2pair_ab_ = local_matrix_kit->vector(dim_ij_ab_);

134

135

Vaa_.restore(si);

136

Vab_.restore(si);

137

Xaa_.restore(si);

138

Xab_.restore(si);

139

Baa_.restore(si);

140

Bab_.restore(si);

141

if (ebc_ == false) {

142

Aaa_.restore(si);

143

Aab_.restore(si);

144

Ac_aa_.restore(si);

145

Ac_ab_.restore(si);

146

T2aa_.restore(si);

147

T2ab_.restore(si);

148

Raa_.restore(si);

149

Rab_.restore(si);

150

}

151

emp2pair_aa_.restore(si);

152

emp2pair_ab_.restore(si);

153

154

int num_tforms;

155

si.get(num_tforms);

156

for(int t=0; t<num_tforms; t++) {

157

std::string tform_name;

158

si.get(tform_name);

159

Ref<TwoBodyMOIntsTransform> tform;

160

tform << SavableState::restore_state(si);

161

tform_map_[tform_name] = tform;

162

}

163

164

int spinadapted; si.get(spinadapted); spinadapted_ = (bool) spinadapted;

165

int evaluated; si.get(evaluated); evaluated_ = (bool) evaluated;

166

si.get(debug_);

167

168

init_tforms_();

169

}

170

171

R12IntEval::~R12IntEval()

172

{

173

r12info_ = 0;

174

dim_ij_aa_ = 0;

175

dim_ij_ab_ = 0;

176

dim_ij_s_ = 0;

177

dim_ij_t_ = 0;

178

dim_ab_aa_ = 0;

179

dim_ab_ab_ = 0;

180

}

181

182

void

183

R12IntEval::save_data_state(StateOut& so)

184

{

185

so.put((int)gbc_);

186

so.put((int)ebc_);

187

so.put((int)abs_method_);

188

so.put((int)stdapprox_);

189

so.put((int)follow_ks_ebcfree_);

190

191

SavableState::save_state(r12info_.pointer(),so);

192

SavableState::save_state(dim_ij_aa_.pointer(),so);

193

SavableState::save_state(dim_ij_ab_.pointer(),so);

194

SavableState::save_state(dim_ij_s_.pointer(),so);

195

SavableState::save_state(dim_ij_t_.pointer(),so);

196

SavableState::save_state(dim_ab_aa_.pointer(),so);

197

SavableState::save_state(dim_ab_ab_.pointer(),so);

198

199

Vaa_.save(so);

200

Vab_.save(so);

201

Xaa_.save(so);

202

Xab_.save(so);

203

Baa_.save(so);

204

Bab_.save(so);

205

if (ebc_ == false) {

206

Aaa_.save(so);

207

Aab_.save(so);

208

Ac_aa_.save(so);

209

Ac_ab_.save(so);

210

T2aa_.save(so);

211

T2ab_.save(so);

212

Raa_.save(so);

213

Rab_.save(so);

214

}

215

emp2pair_aa_.save(so);

216

emp2pair_ab_.save(so);

217

218

int num_tforms = tform_map_.size();

219

so.put(num_tforms);

220

TformMap::iterator first_tform = tform_map_.begin();

221

TformMap::iterator last_tform = tform_map_.end();

222

for(TformMap::iterator t=first_tform; t!=last_tform; t++) {

223

so.put((*t).first);

224

SavableState::save_state((*t).second.pointer(),so);

225

}

226

227

so.put((int)spinadapted_);

228

so.put((int)evaluated_);

229

so.put(debug_);

230

}

231

232

void

233

R12IntEval::obsolete()

234

{

235

evaluated_ = false;

236

237

// make all transforms obsolete

238

TformMap::iterator first_tform = tform_map_.begin();

239

TformMap::iterator last_tform = tform_map_.end();

240

for(TformMap::iterator t=first_tform; t!=last_tform; t++) {

241

(*t).second->obsolete();

242

}

243

244

init_intermeds_();

245

}

246

247

void R12IntEval::include_mp1(bool include_mp1) { include_mp1_ = include_mp1; };

248

void R12IntEval::set_debug(int debug) { if (debug >= 0) { debug_ = debug; r12info_->set_debug_level(debug_); }};

249

void R12IntEval::set_dynamic(bool dynamic) { r12info_->set_dynamic(dynamic); };

250

void R12IntEval::set_print_percent(double pp) { r12info_->set_print_percent(pp); };

251

void R12IntEval::set_memory(size_t nbytes) { r12info_->set_memory(nbytes); };

252

253

Ref<R12IntEvalInfo> R12IntEval::r12info() const { return r12info_; };

254

RefSCDimension R12IntEval::dim_oo_aa() const { return dim_ij_aa_; };

255

RefSCDimension R12IntEval::dim_oo_ab() const { return dim_ij_ab_; };

256

RefSCDimension R12IntEval::dim_oo_s() const { return dim_ij_s_; };

257

RefSCDimension R12IntEval::dim_oo_t() const { return dim_ij_t_; };

258

RefSCDimension R12IntEval::dim_vv_aa() const { return dim_ab_aa_; };

259

RefSCDimension R12IntEval::dim_vv_ab() const { return dim_ab_ab_; };

260

261

RefSCMatrix R12IntEval::V_aa()

262

{

263

compute();

264

return Vaa_;

265

}

266

267

RefSCMatrix R12IntEval::X_aa()

268

{

269

compute();

270

return Xaa_;

271

}

272

273

RefSymmSCMatrix R12IntEval::B_aa()

274

{

275

compute();

276

277

// Extract lower triangle of the matrix

278

Ref<SCMatrixKit> kit = Baa_.kit();

279

RefSymmSCMatrix Baa = kit->symmmatrix(Baa_.rowdim());

280

int naa = Baa_.nrow();

281

double* baa = new double[naa*naa];

282

Baa_.convert(baa);

283

const double* baa_ptr = baa;

284

for(int i=0; i<naa; i++, baa_ptr += i)

285

for(int j=i; j<naa; j++, baa_ptr++)

286

Baa.set_element(i,j,*baa_ptr);

287

delete[] baa;

288

289

return Baa;

290

}

291

292

RefSCMatrix R12IntEval::A_aa()

293

{

294

if (ebc_ == false)

295

compute();

296

return Aaa_;

297

}

298

299

RefSCMatrix R12IntEval::Ac_aa()

300

{

301

if (ebc_ == false && follow_ks_ebcfree_) {

302

compute();

303

return Ac_aa_;

304

}

305

else

306

throw ProgrammingError("R12IntEval::Ac_aa() called although the object initialized with follow_ks_ebcfree set to false",__FILE__,__LINE__);

307

}

308

309

RefSCMatrix R12IntEval::T2_aa()

310

{

311

if (ebc_ == false)

312

compute();

313

return T2aa_;

314

}

315

316

RefSCMatrix R12IntEval::V_ab()

317

{

318

compute();

319

return Vab_;

320

}

321

322

RefSCMatrix R12IntEval::X_ab()

323

{

324

compute();

325

return Xab_;

326

}

327

328

RefSymmSCMatrix R12IntEval::B_ab()

329

{

330

compute();

331

332

// Extract lower triangle of the matrix

333

Ref<SCMatrixKit> kit = Bab_.kit();

334

RefSymmSCMatrix Bab = kit->symmmatrix(Bab_.rowdim());

335

int nab = Bab_.nrow();

336

double* bab = new double[nab*nab];

337

Bab_.convert(bab);

338

const double* bab_ptr = bab;

339

for(int i=0; i<nab; i++, bab_ptr += i)

340

for(int j=i; j<nab; j++, bab_ptr++)

341

Bab.set_element(i,j,*bab_ptr);

342

delete[] bab;

343

344

return Bab;

345

}

346

347

RefSCMatrix R12IntEval::A_ab()

348

{

349

if (ebc_ == false)

350

compute();

351

return Aab_;

352

}

353

354

RefSCMatrix R12IntEval::Ac_ab()

355

{

356

if (ebc_ == false && follow_ks_ebcfree_) {

357

compute();

358

return Ac_ab_;

359

}

360

else

361

throw ProgrammingError("R12IntEval::Ac_ab() called although the object initialized with follow_ks_ebcfree set to false",__FILE__,__LINE__);

362

}

363

364

RefSCMatrix R12IntEval::T2_ab()

365

{

366

if (ebc_ == false)

367

compute();

368

return T2ab_;

369

}

370

371

RefSCVector R12IntEval::emp2_aa()

372

{

373

compute();

374

return emp2pair_aa_;

375

}

376

377

RefSCVector R12IntEval::emp2_ab()

378

{

379

compute();

380

return emp2pair_ab_;

381

}

382

383

RefDiagSCMatrix R12IntEval::evals() const { return r12info_->obs_space()->evals(); };

384

385

void

386

R12IntEval::checkpoint_() const

387

{

388

int me = r12info_->msg()->me();

389

Wavefunction* wfn = r12info_->wfn();

390

391

if (me == 0 && wfn->if_to_checkpoint()) {

392

StateOutBin stateout(wfn->checkpoint_file());

393

SavableState::save_state(wfn,stateout);

394

ExEnv::out0() << indent << "Checkpointed the wave function" << endl;

395

}

396

}

397

398

void

399

R12IntEval::init_tforms_()

400

{

401

Ref<MOIntsTransformFactory> tfactory = r12info_->tfactory();

402

tfactory->set_ints_method((MOIntsTransformFactory::StoreMethod)r12info_->ints_method());

403

404

const std::string ipjq_name = "(ip|jq)";

405

Ref<TwoBodyMOIntsTransform> ipjq_tform = tform_map_[ipjq_name];

406

if (ipjq_tform.null()) {

407

tfactory->set_spaces(r12info_->act_occ_space(),r12info_->obs_space(),

408

r12info_->act_occ_space(),r12info_->obs_space());

409

ipjq_tform = tfactory->twobody_transform_13(ipjq_name);

410

tform_map_[ipjq_name] = ipjq_tform;

411

tform_map_[ipjq_name]->set_num_te_types(3);

412

}

413

414

const std::string iajb_name = "(ia|jb)";

415

Ref<TwoBodyMOIntsTransform> iajb_tform = tform_map_[iajb_name];

416

if (iajb_tform.null()) {

417

tfactory->set_spaces(r12info_->act_occ_space(),r12info_->act_vir_space(),

418

r12info_->act_occ_space(),r12info_->act_vir_space());

419

iajb_tform = tfactory->twobody_transform_13(iajb_name);

420

tform_map_[iajb_name] = iajb_tform;

421

tform_map_[iajb_name]->set_num_te_types(3);

422

}

423

424

const std::string imja_name = "(im|ja)";

425

Ref<TwoBodyMOIntsTransform> imja_tform = tform_map_[imja_name];

426

if (imja_tform.null()) {

427

tfactory->set_spaces(r12info_->act_occ_space(),r12info_->occ_space(),

428

r12info_->act_occ_space(),r12info_->act_vir_space());

429

imja_tform = tfactory->twobody_transform_13(imja_name);

430

tform_map_[imja_name] = imja_tform;

431

tform_map_[imja_name]->set_num_te_types(4);

432

}

433

434

const std::string imjn_name = "(im|jn)";

435

Ref<TwoBodyMOIntsTransform> imjn_tform = tform_map_[imjn_name];

436

if (imjn_tform.null()) {

437

tfactory->set_spaces(r12info_->act_occ_space(),r12info_->occ_space(),

438

r12info_->act_occ_space(),r12info_->occ_space());

439

imjn_tform = tfactory->twobody_transform_13(imjn_name);

440

tform_map_[imjn_name] = imjn_tform;

441

tform_map_[imjn_name]->set_num_te_types(3);

442

}

443

444

const std::string imjy_name = "(im|jy)";

445

Ref<TwoBodyMOIntsTransform> imjy_tform = tform_map_[imjy_name];

446

if (imjy_tform.null()) {

447

tfactory->set_spaces(r12info_->act_occ_space(),r12info_->occ_space(),

448

r12info_->act_occ_space(),r12info_->ribs_space());

449

imjy_tform = tfactory->twobody_transform_13(imjy_name);

450

tform_map_[imjy_name] = imjy_tform;

451

tform_map_[imjy_name]->set_num_te_types(4);

452

}

453

454

iajb_tform = tform_map_[iajb_name];

455

imjn_tform = tform_map_[imjn_name];

456

ipjq_tform = tform_map_[ipjq_name];

457

}

458

459

Ref<TwoBodyMOIntsTransform>

460

R12IntEval::get_tform_(const std::string& tform_name)

461

{

462

TformMap::const_iterator tform_iter = tform_map_.find(tform_name);

463

TwoBodyMOIntsTransform* tform = (*tform_iter).second.pointer();

464

if (tform == NULL) {

465

std::string errmsg = "R12IntEval::get_tform_() -- transform " + tform_name + " is not known";

466

throw std::runtime_error(errmsg.c_str());

467

}

468

tform->compute();

469

470

return tform;

471

}

472

473

void

474

R12IntEval::init_intermeds_()

475

{

476

if (r12info_->msg()->me() == 0) {

477

#if NOT_INCLUDE_DIAGONAL_VXB_CONTIBUTIONS

478

Vaa_->assign(0.0);

479

Vab_->assign(0.0);

480

Baa_->assign(0.0);

481

Bab_->assign(0.0);

482

#else

483

Vaa_->unit();

484

Vab_->unit();

485

Baa_->unit();

486

Bab_->unit();

487

#endif

488

}

489

else {

490

Vaa_.assign(0.0);

491

Vab_.assign(0.0);

492

Baa_.assign(0.0);

493

Bab_.assign(0.0);

494

}

495

if (ebc_ == false) {

496

Aaa_.assign(0.0);

497

Aab_.assign(0.0);

498

if (follow_ks_ebcfree_) {

499

Ac_aa_.assign(0.0);

500

Ac_ab_.assign(0.0);

501

}

502

T2aa_.assign(0.0);

503

T2ab_.assign(0.0);

504

Raa_.assign(0.0);

505

Rab_.assign(0.0);

506

}

507

508

Xaa_.assign(0.0);

509

Xab_.assign(0.0);

510

//r2_contrib_to_X_orig_();

511

#if !NOT_INCLUDE_DIAGONAL_VXB_CONTIBUTIONS

512

r2_contrib_to_X_new_();

513

#endif

514

515

emp2pair_aa_.assign(0.0);

516

emp2pair_ab_.assign(0.0);

517

}

518

519

/// Compute <space1 space1|r_{12}^2|space1 space2>

520

RefSCMatrix

521

R12IntEval::compute_r2_(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2)

522

{

523

/*-----------------------------------------------------

524

Compute overlap, dipole, quadrupole moment integrals

525

-----------------------------------------------------*/

526

RefSCMatrix S_11, MX_11, MY_11, MZ_11, MXX_11, MYY_11, MZZ_11;

527

r12info_->compute_multipole_ints(space1, space1, MX_11, MY_11, MZ_11, MXX_11, MYY_11, MZZ_11);

528

r12info_->compute_overlap_ints(space1, space1, S_11);

529

530

RefSCMatrix S_12, MX_12, MY_12, MZ_12, MXX_12, MYY_12, MZZ_12;

531

if (space1 == space2) {

532

S_12 = S_11;

533

MX_12 = MX_11;

534

MY_12 = MY_11;

535

MZ_12 = MZ_11;

536

MXX_12 = MXX_11;

537

MYY_12 = MYY_11;

538

MZZ_12 = MZZ_11;

539

}

540

else {

541

r12info_->compute_multipole_ints(space1, space2, MX_12, MY_12, MZ_12, MXX_12, MYY_12, MZZ_12);

542

r12info_->compute_overlap_ints(space1, space2, S_12);

543

}

544

if (debug_)

545

ExEnv::out0() << indent << "Computed overlap and multipole moment integrals" << endl;

546

547

const int nproc = r12info_->msg()->n();

548

const int me = r12info_->msg()->me();

549

550

const int n1 = space1->rank();

551

const int n2 = space2->rank();

552

const int n12 = n1*n2;

553

const int n1112 = n1*n1*n12;

554

double* r2_array = new double[n1112];

555

memset(r2_array,0,n1112*sizeof(double));

556

557

int ij = 0;

558

double* ijkl_ptr = r2_array;

559

for(int i=0; i<n1; i++)

560

for(int j=0; j<n1; j++, ij++) {

561

562

int ij_proc = ij%nproc;

563

if (ij_proc != me) {

564

ijkl_ptr += n12;

565

continue;

566

}

567

568

int kl=0;

569

for(int k=0; k<n1; k++)

570

for(int l=0; l<n2; l++, kl++, ijkl_ptr++) {

571

572

double r1r1_ik = -1.0*(MXX_11->get_element(i,k) + MYY_11->get_element(i,k) + MZZ_11->get_element(i,k));

573

double r1r1_jl = -1.0*(MXX_12->get_element(j,l) + MYY_12->get_element(j,l) + MZZ_12->get_element(j,l));

574

double r1r2_ijkl = MX_11->get_element(i,k)*MX_12->get_element(j,l) +

575

MY_11->get_element(i,k)*MY_12->get_element(j,l) +

576

MZ_11->get_element(i,k)*MZ_12->get_element(j,l);

577

double S_ik = S_11.get_element(i,k);

578

double S_jl = S_12.get_element(j,l);

579

580

double R2_ijkl = r1r1_ik * S_jl + r1r1_jl * S_ik - 2.0*r1r2_ijkl;

581

*ijkl_ptr = R2_ijkl;

582

}

583

}

584

585

r12info_->msg()->sum(r2_array,n1112);

586

587

MOPairIterFactory pair_factory;

588

RefSCDimension dim_ij = pair_factory.scdim_ab(space1,space1);

589

RefSCDimension dim_kl = pair_factory.scdim_ab(space1,space2);

590

591

Ref<LocalSCMatrixKit> local_matrix_kit = new LocalSCMatrixKit();

592

RefSCMatrix R2 = local_matrix_kit->matrix(dim_ij, dim_kl);

593

R2.assign(r2_array);

594

delete[] r2_array;

595

596

return R2;

597

}

598

599

600

void

601

R12IntEval::r2_contrib_to_X_orig_()

602

{

603

/*---------------------------------------------------------------

604

Compute dipole and quadrupole moment integrals in act MO basis

605

---------------------------------------------------------------*/

606

RefSCMatrix MX, MY, MZ, MXX, MYY, MZZ;

607

r12info_->compute_multipole_ints(r12info_->act_occ_space(),r12info_->act_occ_space(),MX,MY,MZ,MXX,MYY,MZZ);

608

if (debug_)

609

ExEnv::out0() << indent << "Computed multipole moment integrals" << endl;

610

611

const int nproc = r12info_->msg()->n();

612

const int me = r12info_->msg()->me();

613

614

SpatialMOPairIter_eq ij_iter(r12info_->act_occ_space());

615

SpatialMOPairIter_eq kl_iter(r12info_->act_occ_space());

616

617

for(kl_iter.start();int(kl_iter);kl_iter.next()) {

618

619

const int kl = kl_iter.ij();

620

int kl_proc = kl%nproc;

621

if (kl_proc != me)

622

continue;

623

const int k = kl_iter.i();

624

const int l = kl_iter.j();

625

const int kl_aa = kl_iter.ij_aa();

626

const int kl_ab = kl_iter.ij_ab();

627

const int lk_ab = kl_iter.ij_ba();

628

629

for(ij_iter.start();int(ij_iter);ij_iter.next()) {

630

631

const int i = ij_iter.i();

632

const int j = ij_iter.j();

633

const int ij_aa = ij_iter.ij_aa();

634

const int ij_ab = ij_iter.ij_ab();

635

const int ji_ab = ij_iter.ij_ba();

636

637

/*----------------------------------

638

Compute (r12)^2 contribution to X

639

----------------------------------*/

640

double r1r1_ik = -1.0*(MXX->get_element(i,k) + MYY->get_element(i,k) + MZZ->get_element(i,k));

641

double r1r1_il = -1.0*(MXX->get_element(i,l) + MYY->get_element(i,l) + MZZ->get_element(i,l));

642

double r1r1_jk = -1.0*(MXX->get_element(j,k) + MYY->get_element(j,k) + MZZ->get_element(j,k));

643

double r1r1_jl = -1.0*(MXX->get_element(j,l) + MYY->get_element(j,l) + MZZ->get_element(j,l));

644

double r1r2_ijkl = MX->get_element(i,k)*MX->get_element(j,l) +

645

MY->get_element(i,k)*MY->get_element(j,l) +

646

MZ->get_element(i,k)*MZ->get_element(j,l);

647

double r1r2_ijlk = MX->get_element(i,l)*MX->get_element(j,k) +

648

MY->get_element(i,l)*MY->get_element(j,k) +

649

MZ->get_element(i,l)*MZ->get_element(j,k);

650

double delta_ik = (i==k ? 1.0 : 0.0);

651

double delta_il = (i==l ? 1.0 : 0.0);

652

double delta_jk = (j==k ? 1.0 : 0.0);

653

double delta_jl = (j==l ? 1.0 : 0.0);

654

655

double Xab_ijkl = r1r1_ik * delta_jl + r1r1_jl * delta_ik - 2.0*r1r2_ijkl;

656

Xab_.accumulate_element(ij_ab,kl_ab,Xab_ijkl);

657

658

if (ij_ab != ji_ab) {

659

double Xab_jikl = r1r1_jk * delta_il + r1r1_il * delta_jk - 2.0*r1r2_ijlk;

660

Xab_.accumulate_element(ji_ab,kl_ab,Xab_jikl);

661

}

662

663

if (kl_ab != lk_ab) {

664

double Xab_ijlk = r1r1_il * delta_jk + r1r1_jk * delta_il - 2.0*r1r2_ijlk;

665

Xab_.accumulate_element(ij_ab,lk_ab,Xab_ijlk);

666

}

667

668

if (ij_ab != ji_ab && kl_ab != lk_ab) {

669

double Xab_jilk = r1r1_ik * delta_jl + r1r1_jl * delta_ik - 2.0*r1r2_ijkl;

670

Xab_.accumulate_element(ji_ab,lk_ab,Xab_jilk);

671

}

672

673

if (ij_aa != -1 && kl_aa != -1) {

674

double Xaa_ijkl = r1r1_ik * delta_jl + r1r1_jl * delta_ik - 2.0*r1r2_ijkl -

675

r1r1_jk * delta_il - r1r1_il * delta_jk + 2.0*r1r2_ijlk;

676

Xaa_.accumulate_element(ij_aa,kl_aa,Xaa_ijkl);

677

}

678

679

}

680

}

681

}

682

683

void

684

R12IntEval::r2_contrib_to_X_new_()

685

{

686

unsigned int me = r12info_->msg()->me();

687

688

// compute r_{12}^2 operator in act.occ.pair/act.occ.pair basis

689

RefSCMatrix R2 = compute_r2_(r12info_->act_occ_space(),r12info_->act_occ_space());

690

691

if (me != 0)

692

return;

693

Xab_.accumulate(R2);

694

695

SpatialMOPairIter_eq ij_iter(r12info_->act_occ_space());

696

SpatialMOPairIter_eq kl_iter(r12info_->act_occ_space());

697

698

for(kl_iter.start();int(kl_iter);kl_iter.next()) {

699

700

const int kl_aa = kl_iter.ij_aa();

701

if (kl_aa == -1)

702

continue;

703

const int kl_ab = kl_iter.ij_ab();

704

705

for(ij_iter.start();int(ij_iter);ij_iter.next()) {

706

707

const int ij_aa = ij_iter.ij_aa();

708

const int ij_ab = ij_iter.ij_ab();

709

const int ji_ab = ij_iter.ij_ba();

710

711

if (ij_aa != -1) {

712

double Xaa_ijkl = R2.get_element(ij_ab,kl_ab) - R2.get_element(ji_ab,kl_ab);

713

Xaa_.accumulate_element(ij_aa,kl_aa,Xaa_ijkl);

714

}

715

716

}

717

}

718

}

719

720

void

721

R12IntEval::form_focc_space_()

722

{

723

// compute the Fock matrix between the complement and all occupieds and

724

// create the new Fock-weighted space

725

if (focc_space_.null()) {

726

Ref<MOIndexSpace> occ_space = r12info_->occ_space();

727

Ref<MOIndexSpace> ribs_space = r12info_->ribs_space();

728

729

RefSCMatrix F_ri_o = fock_(occ_space,ribs_space,occ_space);

730

if (debug_ > 1)

731

F_ri_o.print("Fock matrix (RI-BS/occ.)");

732

focc_space_ = new MOIndexSpace("Fock-weighted occupied MOs sorted by energy",

733

occ_space, ribs_space->coefs()*F_ri_o, ribs_space->basis());

734

}

735

}

736

737

void

738

R12IntEval::form_factocc_space_()

739

{

740

// compute the Fock matrix between the complement and active occupieds and

741

// create the new Fock-weighted space

742

if (factocc_space_.null()) {

743

Ref<MOIndexSpace> occ_space = r12info_->occ_space();

744

Ref<MOIndexSpace> act_occ_space = r12info_->act_occ_space();

745

Ref<MOIndexSpace> ribs_space = r12info_->ribs_space();

746

747

RefSCMatrix F_ri_ao = fock_(occ_space,ribs_space,act_occ_space);

748

if (debug_ > 1)

749

F_ri_ao.print("Fock matrix (RI-BS/act.occ.)");

750

factocc_space_ = new MOIndexSpace("Fock-weighted active occupied MOs sorted by energy",

751

act_occ_space, ribs_space->coefs()*F_ri_ao, ribs_space->basis());

752

}

753

}

754

755

void

756

R12IntEval::form_canonvir_space_()

757

{

758

// Create a complement space to all occupieds

759

// Fock operator is diagonal in this space

760

if (canonvir_space_.null()) {

761

762

if (r12info_->basis_vir()->equiv(r12info_->basis())) {

763

canonvir_space_ = r12info_->vir_space();

764

return;

765

}

766

767

const Ref<MOIndexSpace> mo_space = r12info_->mo_space();

768

Ref<MOIndexSpace> vir_space = r12info_->vir_space_symblk();

769

RefSCMatrix F_vir = fock_(r12info_->occ_space(),vir_space,vir_space);

770

771

int nrow = vir_space->rank();

772

double* F_full = new double[nrow*nrow];

773

double* F_lowtri = new double [nrow*(nrow+1)/2];

774

F_vir->convert(F_full);

775

int ij = 0;

776

for(int row=0; row<nrow; row++) {

777

int rc = row*nrow;

778

for(int col=0; col<=row; col++, rc++, ij++) {

779

F_lowtri[ij] = F_full[rc];

780

}

781

}

782

RefSymmSCMatrix F_vir_lt(F_vir.rowdim(),F_vir->kit());

783

F_vir_lt->assign(F_lowtri);

784

F_vir = 0;

785

delete[] F_full;

786

delete[] F_lowtri;

787

788

Ref<MOIndexSpace> canonvir_space_symblk = new MOIndexSpace("Virt. MOs symmetry-blocked",

789

vir_space, vir_space->coefs()*F_vir_lt.eigvecs(),

790

vir_space->basis());

791

RefDiagSCMatrix F_vir_evals = F_vir_lt.eigvals();

792

canonvir_space_ = new MOIndexSpace("Virt. MOs sorted by energy",

793

canonvir_space_symblk->coefs(),

794

canonvir_space_symblk->basis(),

795

canonvir_space_symblk->integral(),

796

F_vir_evals, 0, 0,

797

MOIndexSpace::energy);

798

}

799

}

800

801

const int

802

R12IntEval::tasks_with_ints_(const Ref<R12IntsAcc> ints_acc, vector<int>& map_to_twi)

803

{

804

int nproc = r12info_->msg()->n();

805

806

// Compute the number of tasks that have full access to the integrals

807

// and split the work among them

808

int nproc_with_ints = 0;

809

for(int proc=0;proc<nproc;proc++)

810

if (ints_acc->has_access(proc)) nproc_with_ints++;

811

812

map_to_twi.resize(nproc);

813

int count = 0;

814

for(int proc=0;proc<nproc;proc++)

815

if (ints_acc->has_access(proc)) {

816

map_to_twi[proc] = count;

817

count++;

818

}

819

else

820

map_to_twi[proc] = -1;

821

822

ExEnv::out0() << indent << "Computing intermediates on " << nproc_with_ints

823

<< " processors" << endl;

824

825

return nproc_with_ints;

826

}

827

828

829

void

830

R12IntEval::compute()

831

{

832

if (evaluated_)

833

return;

834

835

if (r12info_->basis_vir()->equiv(r12info_->basis())) {

836

obs_contrib_to_VXB_gebc_vbseqobs_();

837

if (debug_ > 1) {

838

Vaa_.print("Alpha-alpha V(OBS) contribution");

839

Vab_.print("Alpha-beta V(OBS) contribution");

840

Xaa_.print("Alpha-alpha X(OBS) contribution");

841

Xab_.print("Alpha-beta X(OBS) contribution");

842

Baa_.print("Alpha-alpha B(OBS) contribution");

843

Bab_.print("Alpha-beta B(OBS) contribution");

844

}

845

if (r12info_->basis() != r12info_->basis_ri())

846

abs1_contrib_to_VXB_gebc_();

847

if (debug_ > 1) {

848

Vaa_.print("Alpha-alpha V(OBS+ABS) contribution");

849

Vab_.print("Alpha-beta V(OBS+ABS) contribution");

850

Xaa_.print("Alpha-alpha X(OBS+ABS) contribution");

851

Xab_.print("Alpha-beta X(OBS+ABS) contribution");

852

Baa_.print("Alpha-alpha B(OBS+ABS) contribution");

853

Bab_.print("Alpha-beta B(OBS+ABS) contribution");

854

}

855

}

856

else {

857

contrib_to_VXB_gebc_vbsneqobs_();

858

compute_dualEmp2_();

859

if (include_mp1_)

860

compute_dualEmp1_();

861

}

862

863

864

#if TEST_FOCK

865

if (!evaluated_) {

866

RefSCMatrix F = fock_(r12info_->occ_space(),r12info_->obs_space(),r12info_->obs_space());

867

F.print("Fock matrix in OBS");

868

r12info_->obs_space()->evals().print("OBS eigenvalues");

869

870

r12info_->ribs_space()->coefs().print("Orthonormal RI-BS");

871

RefSCMatrix S_ri;

872

r12info_->compute_overlap_ints(r12info_->ribs_space(),r12info_->ribs_space(),S_ri);

873

S_ri.print("Overlap in RI-BS");

874

RefSCMatrix F_ri = fock_(r12info_->occ_space(),r12info_->ribs_space(),r12info_->ribs_space());

875

F_ri.print("Fock matrix in RI-BS");

876

RefSymmSCMatrix F_ri_symm = F_ri.kit()->symmmatrix(F_ri.rowdim());

877

int nrow = F_ri.rowdim().n();

878

for(int r=0; r<nrow; r++)

879

for(int c=0; c<nrow; c++)

880

F_ri_symm.set_element(r,c,F_ri.get_element(r,c));

881

F_ri_symm.eigvals().print("Eigenvalues of the Fock matrix (RI-BS)");

882

883

RefSCMatrix F_obs_ri = fock_(r12info_->occ_space(),r12info_->obs_space(),r12info_->ribs_space());

884

F_obs_ri.print("Fock matrix in OBS/RI-BS");

885

}

886

#endif

887

888

if (!ebc_) {

889

// These functions assume that virtuals are expanded in the same basis

890

// as the occupied orbitals

891

if (!r12info_->basis_vir()->equiv(r12info_->basis()))

892

throw std::runtime_error("R12IntEval::compute() -- ebc=false is only supported when basis_vir == basis");

893

894

compute_A_simple_();

895

Aaa_.scale(2.0);

896

Aab_.scale(2.0);

897

if (follow_ks_ebcfree_) {

898

compute_A_via_commutator_();

899

Ac_aa_.scale(2.0);

900

Ac_ab_.scale(2.0);

901

}

902

compute_T2_();

903

AT2_contrib_to_V_();

904

compute_R_();

905

AR_contrib_to_B_();

906

}

907

908

if (!gbc_) {

909

// These functions assume that virtuals are expanded in the same basis

910

// as the occupied orbitals

911

if (!r12info_->basis_vir()->equiv(r12info_->basis()))

912

throw std::runtime_error("R12IntEval::compute() -- gbc=false is only supported when basis_vir == basis");

913

914

compute_B_gbc_1_();

915

if (debug_ > 1) {

916

Baa_.print("Alpha-alpha B(OBS+ABS+GBC1) contribution");

917

Bab_.print("Alpha-beta B(OBS+ABS+GBC1) contribution");

918

}

919

compute_B_gbc_2_();

920

if (debug_ > 1) {

921

Baa_.print("Alpha-alpha B(OBS+ABS+GBC1+GBC2) contribution");

922

Bab_.print("Alpha-beta B(OBS+ABS+GBC1+GBC2) contribution");

923

}

924

}

925

926

// Distribute the final intermediates to every node

927

globally_sum_intermeds_(true);

928

929

evaluated_ = true;

930

}

931

932

void

933

R12IntEval::globally_sum_scmatrix_(RefSCMatrix& A, bool to_all_tasks, bool to_average)

934

{

935

Ref<MessageGrp> msg = r12info_->msg();

936

unsigned int ntasks = msg->n();

937

// If there's only one task then there's nothing to do

938

if (ntasks == 1)

939

return;

940

941

const int nelem = A.ncol() * A.nrow();

942

double *A_array = new double[nelem];

943

A.convert(A_array);

944

if (to_all_tasks)

945

msg->sum(A_array,nelem,0,-1);

946

else

947

msg->sum(A_array,nelem,0,0);

948

A.assign(A_array);

949

if (to_average)

950

A.scale(1.0/(double)ntasks);

951

if (!to_all_tasks && msg->me() != 0)

952

A.assign(0.0);

953

954

delete[] A_array;

955

}

956

957

void

958

R12IntEval::globally_sum_scvector_(RefSCVector& A, bool to_all_tasks, bool to_average)

959

{

960

Ref<MessageGrp> msg = r12info_->msg();

961

unsigned int ntasks = msg->n();

962

// If there's only one task then there's nothing to do

963

if (ntasks == 1)

964

return;

965

966

const int nelem = A.dim().n();

967

double *A_array = new double[nelem];

968

A.convert(A_array);

969

if (to_all_tasks)

970

msg->sum(A_array,nelem,0,-1);

971

else

972

msg->sum(A_array,nelem,0,0);

973

A.assign(A_array);

974

if (to_average)

975

A.scale(1.0/(double)ntasks);

976

if (!to_all_tasks && msg->me() != 0)

977

A.assign(0.0);

978

979

delete[] A_array;

980

}

981

982

void

983

R12IntEval::globally_sum_intermeds_(bool to_all_tasks)

984

{

985

globally_sum_scmatrix_(Vaa_,to_all_tasks);

986

globally_sum_scmatrix_(Vab_,to_all_tasks);

987

988

globally_sum_scmatrix_(Xaa_,to_all_tasks);

989

globally_sum_scmatrix_(Xab_,to_all_tasks);

990

991

globally_sum_scmatrix_(Baa_,to_all_tasks);

992

globally_sum_scmatrix_(Bab_,to_all_tasks);

993

994

if (ebc_ == false) {

995

globally_sum_scmatrix_(Aaa_,to_all_tasks);

996

globally_sum_scmatrix_(Aab_,to_all_tasks);

997

998

if (follow_ks_ebcfree_) {

999

globally_sum_scmatrix_(Ac_aa_,to_all_tasks);

1000

globally_sum_scmatrix_(Ac_ab_,to_all_tasks);

1001

}

1002

1003

globally_sum_scmatrix_(T2aa_,to_all_tasks);

1004

globally_sum_scmatrix_(T2ab_,to_all_tasks);

1005

1006

globally_sum_scmatrix_(Raa_,to_all_tasks);

1007

globally_sum_scmatrix_(Rab_,to_all_tasks);

1008

}

1009

1010

globally_sum_scvector_(emp2pair_aa_,to_all_tasks);

1011

globally_sum_scvector_(emp2pair_ab_,to_all_tasks);

1012

1013

if (debug_) {

1014

ExEnv::out0() << indent << "Collected contributions to the intermediates from all tasks";

1015

if (to_all_tasks)

1016

ExEnv::out0() << " and distributed to every task" << endl;

1017

else

1018

ExEnv::out0() << " on task 0" << endl;

1019

}

1020

}

1021

1022

/////////////////////////////////////////////////////////////////////////////

1023

1024

// Local Variables:

1025

// mode: c++

1026

// c-file-style: "CLJ-CONDENSED"

1027

// End:

Older »