3
\brief Enter brief description of file here
6
* Revision 1.6 2004/05/03 04:32:40 crawdad
7
* Major mods based on merge with stable psi-3-2-1 release. Note that this
8
* version has not been fully tested and some scf-optn test cases do not run
9
* correctly beccause of changes in mid-March 2004 to optking.
12
/* Revision 1.5.8.1 2004/04/10 19:41:32 crawdad
13
/* Fixed the DIIS code for UHF cases. The new version uses the Pulay scheme of
14
/* building the error vector in the AO basis as FDS-SDF, followed by xformation
15
/* into the orthogonal AO basis. This code converges faster for test cases
16
/* like cc8, but fails for linearly dependent basis sets for unknown reasons.
19
/* Revision 1.5 2002/04/03 02:06:01 janssen
20
/* Finish changes to use new include paths for libraries.
22
/* Revision 1.4 2000/12/05 19:40:03 sbrown
23
/* Added Unrestricted Kohn-Sham DFT.
25
/* Revision 1.3 2000/06/22 22:15:00 evaleev
26
/* Modifications for KS DFT. Reading in XC Fock matrices and XC energy in formg_direct need to be uncommented (at present those are not produced by CINTS yet).
28
/* Revision 1.2 2000/06/02 13:32:15 kenny
31
/* Added dynamic integral accuracy cutoffs for direct scf. Added a few global
32
/* variables. Added keyword 'dyn_acc'; true--use dynamic cutoffs. Use of
33
/* 'dconv' and 'delta' to keep track of density convergence somewhat awkward,
34
/* but avoids problems when accuracy is switched and we have to wipe out density
35
/* matrices. Also added error message and exit if direct rohf singlet is
36
/* attempted since it doesn't work.
39
/* Revision 1.1.1.1 2000/02/04 22:52:33 evaleev
40
/* Started PSI 3 repository
42
/* Revision 1.2 1999/11/17 19:40:46 evaleev
43
/* Made all the adjustments necessary to have direct UHF working. Still doesn't work though..
45
/* Revision 1.1 1999/11/02 23:55:56 sbrown
46
/* Shawn Brown - (11/2/99) Modified to the code in a few major ways.
48
/* 1. Added the capability to do UHF. All of the features available with the
49
/* other refrences have been added for UHF.
51
/* 2. For UHF, I had to alter the structure of file30. (See cleanup.c for a
52
/* map) This entailed adding a pointer array right after the header in the SCF
53
/* section of file30 that pointed to all of the data for the SCF caclulation.
54
/* Functions were added to libfile30 to account for this and they are
55
/* incorporated in this code.
57
/* 3. Updated and fixed all of the problems associated with my previous
58
/* guessing code. The code no longer uses OPENTYPE to specify the type of
59
/* occupation. The keword REFERENCE and MULTP can now be used to indicate any
60
/* type of calculation. (e.g. ROHF with MULTP of 1 is an open shell singlet
61
/* ROHF calculation) This code was moved to occ_fun.c. The code can also
62
/* guess at any multplicity in a highspin case, provided enough electrons.
64
/* Revision 1.1.1.1 1999/04/12 16:59:25 evaleev
65
/* Added a version of CSCF that can work with CINTS.
68
* Revision 1.1 1991/06/15 20:22:21 seidl
73
#include <libpsio/psio.h>
77
namespace psi { namespace cscf {
81
int i,j,k,l,ij,n,m,jj,kk;
91
for(m = 0; m < 2; m++){
92
for (l=0; l < num_ir ; l++) {
94
if (n=scf_info[l].num_so) {
95
ndocc = sp->scf_spin[l].noccup;
96
for (i=ij=0; i < n ; i++ ) {
97
/*--------------------------------------
101
-------------------------------------*/
103
for (j=0; j <i; j++,ij++) {
105
for (k=0; k < ndocc ; k++)
106
ptemp += 2.0*sp->scf_spin[l].cmat[i][k]
107
*sp->scf_spin[l].cmat[j][k];
108
sp->scf_spin[l].dpmat[ij] = ptemp - sp->scf_spin[l].pmat[ij];
109
sp->scf_spin[l].pmat[ij] = ptemp;
111
/*----------------------------------
115
--------------------------------*/
117
for (k=0; k < ndocc ; k++) {
118
ctmp=sp->scf_spin[l].cmat[i][k];
121
sp->scf_spin[l].dpmat[ij] = ptemp - sp->scf_spin[l].pmat[ij];
122
sp->scf_spin[l].pmat[ij] = ptemp;
127
"\nSpin case %d density matrix for irrep %s",m,scf_info[l].irrep_label);
128
print_array(spin_info[m].scf_spin[l].pmat,
129
scf_info[l].num_so,outfile);
135
for(l=0;l < num_ir; l++){
136
if(nn=scf_info[l].num_so) {
137
for(ij=0;ij<ioff[nn];ij++) {
138
ptemp = spin_info[0].scf_spin[l].pmat[ij] +
139
spin_info[1].scf_spin[l].pmat[ij];
140
scf_info[l].dpmat[ij] = ptemp - scf_info[l].pmat[ij];
141
scf_info[l].pmat[ij] = ptemp;
146
/*-----------------------
147
Handle direct SCF here
148
-----------------------*/
150
/*decide what accuracy to request for direct_scf*/
152
if((iter<30)&&(tight_ints==0)&&(delta>1.0E-5)) {
155
if((tight_ints==0)&&(delta<=1.0E-5)){
156
fprintf(outfile," Switching to full integral accuracy\n");
163
psio_open(itapDSCF, PSIO_OPEN_NEW);
164
psio_write_entry(itapDSCF, "Integrals cutoff", (char *) &eri_cutoff, sizeof(double));
166
ntri = nbasis*(nbasis+1)/2;
167
dmat = init_array(ntri);
169
/*--- Get full dpmata ---*/
170
for(i=0;i<num_ir;i++) {
171
max = scf_info[i].num_so;
172
off = scf_info[i].ideg;
178
dmat[ioff[jj]+kk] = spin_info[0].scf_spin[i].pmat[ioff[j]+k];
179
spin_info[0].scf_spin[i].dpmat[ioff[j]+k] = 0.0;
182
dmat[ioff[jj]+kk] = spin_info[0].scf_spin[i].dpmat[ioff[j]+k];
187
psio_write_entry(itapDSCF, "Difference Alpha Density", (char *) dmat, sizeof(double)*ntri);
189
/*--- Get full dpmatb ---*/
190
for(i=0;i<num_ir;i++) {
191
max = scf_info[i].num_so;
192
off = scf_info[i].ideg;
198
dmat[ioff[jj]+kk] = spin_info[1].scf_spin[i].pmat[ioff[j]+k];
199
spin_info[1].scf_spin[i].dpmat[ioff[j]+k] = 0.0;
202
dmat[ioff[jj]+kk] = spin_info[1].scf_spin[i].dpmat[ioff[j]+k];
207
psio_write_entry(itapDSCF, "Difference Beta Density", (char *) dmat, sizeof(double)*ntri);
210
/* ---- Get Occupied Eigenvector matrix for DFT ----*/
212
/*------ Get the Alpha Occupied Eigenvector --------*/
214
cmat = block_matrix(nbfso,a_elec);
215
for(i=l=0;i<num_ir;i++){
216
max = spin_info[0].scf_spin[i].nclosed;
217
off = scf_info[i].ideg;
219
for(k=0;k<scf_info[i].num_mo;k++) {
221
cmat[kk][l] = spin_info[0].scf_spin[i].cmat[k][j];
226
/*fprintf(outfile,"\na_elec = %d",a_elec);
227
fprintf(outfile,"\nAlpha Occupied Eigenvector from CSCF");
228
print_mat(cmat,nbfso,a_elec,outfile);*/
230
psio_write_entry(itapDSCF, "Number of MOs",
231
(char *) &(nmo),sizeof(int));
233
psio_write_entry(itapDSCF, "Number of Alpha DOCC",
234
(char *) &(a_elec),sizeof(int));
236
psio_write_entry(itapDSCF, "Alpha Occupied SCF Eigenvector",
237
(char *) &(cmat[0][0]),sizeof(double)*ntri);
238
/*fprintf(outfile,"\na_elec = %d",a_elec);*/
240
/*------ Get the Alpha Occupied Eigenvector --------*/
242
cmat = block_matrix(nbfso,b_elec);
243
for(i=l=0;i<num_ir;i++){
244
max = spin_info[1].scf_spin[i].nclosed;
245
off = scf_info[i].ideg;
247
for(k=0;k<scf_info[i].num_mo;k++) {
249
cmat[kk][l] = spin_info[1].scf_spin[i].cmat[k][j];
254
/*fprintf(outfile,"\nBeta Occupied Eigenvector from CSCF");
255
print_mat(cmat,nbfso,b_elec,outfile);*/
257
psio_write_entry(itapDSCF, "Number of Beta DOCC",
258
(char *) &(b_elec),sizeof(int));
259
psio_write_entry(itapDSCF, "Beta Occupied SCF Eigenvector",
260
(char *) &(cmat[0][0]),sizeof(double)*ntri);
263
/*--- Get full dpmata ---*/
264
for(i=0;i<num_ir;i++) {
265
max = scf_info[i].num_so;
266
off = scf_info[i].ideg;
271
dmat[ioff[jj]+kk] = spin_info[0].scf_spin[i].pmat[ioff[j]+k];
275
psio_write_entry(itapDSCF, "Total Alpha Density", (char *) dmat, sizeof(double)*ntri);
277
/*--- Get full dpmatb ---*/
278
for(i=0;i<num_ir;i++) {
279
max = scf_info[i].num_so;
280
off = scf_info[i].ideg;
285
dmat[ioff[jj]+kk] = spin_info[1].scf_spin[i].pmat[ioff[j]+k];
289
psio_write_entry(itapDSCF, "Total Beta Density", (char *) dmat, sizeof(double)*ntri);
292
psio_close(itapDSCF, 1);
298
}} // namespace psi::cscf