2
** CCLAMBDA: Program to calculate the coupled-cluster lambda vector.
8
#include <libipv1/ip_lib.h>
9
#include <libpsio/psio.h>
10
#include <libciomr/libciomr.h>
11
#include <libdpd/dpd.h>
15
/* Function prototypes */
16
void init_io(int argc, char *argv[]);
18
void get_moinfo(void);
19
void get_params(void);
21
void init_amps(int L_irr, int state__index);
22
double pseudoenergy(int L_irr);
24
void G_build(int L_irr);
25
void L1_build(int L_irr, int root_L_irr);
26
void L2_build(int L_irr, int root_L_irr);
27
void sort_amps(int L_irr);
28
void Lsave(int L_irr);
29
void Lnorm(int L_irr, int root_L_irr);
32
int converged(int L_irr);
33
int **cacheprep_rhf(int level, int *cachefiles);
34
int **cacheprep_uhf(int level, int *cachefiles);
35
void cachedone_rhf(int **cachelist);
36
void cachedone_uhf(int **cachelist);
37
void denom(int L_irr, int e_index);
38
void overlap(int L_irr);
39
void Lsave_index(int L_irr, int root_L_irr);
40
void L_to_LAMPS(int L_irr);
41
extern void check_ortho(void);
42
void L_zero(int L_irr);
44
int main(int argc, char *argv[])
46
int done=0, i, L_irr, root_L_irr;
47
int **cachelist, *cachefiles;
50
init_io(argc, argv); /* parses command-line arguments */
56
cachefiles = init_int_array(PSIO_MAXUNIT);
58
if(params.ref == 0 || params.ref == 1) { /** RHF or ROHF **/
60
cachelist = cacheprep_rhf(params.cachelev, cachefiles);
62
dpd_init(0, moinfo.nirreps, params.memory, 0, cachefiles, cachelist, NULL,
63
2, moinfo.occpi, moinfo.occ_sym, moinfo.virtpi, moinfo.vir_sym);
65
if(params.aobasis) { /* Set up new DPD for AO-basis algorithm */
66
dpd_init(1, moinfo.nirreps, params.memory, 0, cachefiles, cachelist, NULL,
67
2, moinfo.occpi, moinfo.orbsym, moinfo.orbspi, moinfo.orbsym);
72
else if(params.ref == 2) { /** UHF **/
74
cachelist = cacheprep_uhf(params.cachelev, cachefiles);
76
dpd_init(0, moinfo.nirreps, params.memory, 0, cachefiles,
77
cachelist, NULL, 4, moinfo.aoccpi, moinfo.aocc_sym, moinfo.avirtpi,
78
moinfo.avir_sym, moinfo.boccpi, moinfo.bocc_sym, moinfo.bvirtpi, moinfo.bvir_sym);
80
if(params.aobasis) { /* Set up new DPD's for AO-basis algorithm */
81
dpd_init(1, moinfo.nirreps, params.memory, 0, cachefiles, cachelist, NULL,
82
4, moinfo.aoccpi, moinfo.aocc_sym, moinfo.orbspi, moinfo.orbsym,
83
moinfo.boccpi, moinfo.bocc_sym, moinfo.orbspi, moinfo.orbsym);
89
if(params.local) local_init();
91
for (L_irr=0; L_irr<moinfo.nirreps; ++L_irr) {
93
psio_close(CC_TMP1,0);
94
psio_close(CC_TMP2,0);
95
psio_close(CC_LAMBDA,0);
99
psio_open(CC_LAMBDA,0);
100
for (root_L_irr=0; root_L_irr < params.Ls_per_irrep[L_irr]; ++root_L_irr) {
101
fprintf(outfile,"\tSymmetry of excited state: %s\n", moinfo.labels[moinfo.sym^L_irr]);
102
fprintf(outfile,"\tSymmetry of right eigenvector: %s\n",moinfo.labels[L_irr]);
104
init_amps(L_irr, root_L_irr);
106
fprintf(outfile, "\n\t Solving Lambda Equations\n");
107
fprintf(outfile, "\t ------------------------\n");
108
fprintf(outfile, "\tIter PseudoEnergy or Norm RMS \n");
109
fprintf(outfile, "\t---- --------------------- --------\n");
111
denom(L_irr, root_L_irr); /* second argument determines E used */
113
moinfo.lcc = pseudoenergy(L_irr);
116
for(moinfo.iter=1 ; moinfo.iter <= params.maxiter; moinfo.iter++) {
120
/* must zero New L before adding RHS */
122
L1_build(L_irr,root_L_irr);
123
L2_build(L_irr,root_L_irr);
125
if(converged(L_irr)) {
126
done = 1; /* Boolean for convergence */
127
Lsave(L_irr); /* copy "New L" to "L" */
128
moinfo.lcc = pseudoenergy(L_irr);
130
if (!params.ground) {
131
Lnorm(L_irr, root_L_irr); /* normalize against R */
132
Lsave_index(L_irr, root_L_irr); /* put Ls in unique location */
137
/* sort_amps(); to be done by later functions */
138
fprintf(outfile, "\n\tIterations converged.\n");
143
diis(moinfo.iter, L_irr);
145
moinfo.lcc = pseudoenergy(L_irr);
148
fprintf(outfile, "\n");
150
fprintf(outfile, "\t ** Lambda not converged to %2.1e ** \n",
156
exit(PSI_RETURN_FAILURE);
160
overlap_LAMPS(L_irr);
165
if(params.local) local_done();
167
if (!params.ground) check_ortho();
171
if(params.ref == 2) cachedone_uhf(cachelist);
172
else cachedone_rhf(cachelist);
177
exit(PSI_RETURN_SUCCESS);
180
/* parse command line arguments */
181
void init_io(int argc, char *argv[])
184
extern char *gprgid();
185
char *lbl, *progid, *argv_unparsed[100];
187
progid = (char *) malloc(strlen(gprgid())+2);
188
sprintf(progid, ":%s",gprgid());
191
for (i=1, num_unparsed=0; i<argc; ++i) {
192
if (!strcmp(argv[i],"--excited")) {
196
argv_unparsed[num_unparsed++] = argv[i];
200
psi_start(num_unparsed, argv_unparsed, 0);
201
ip_cwk_add(":INPUT");
205
ip_cwk_add(":CCEOM");
210
for(i=CC_MIN; i <= CC_MAX; i++) psio_open(i,1);
221
fprintf(outfile, "\n");
222
fprintf(outfile, "\t\t\t**************************\n");
223
fprintf(outfile, "\t\t\t* CCLAMBDA *\n");
224
fprintf(outfile, "\t\t\t**************************\n");
225
fprintf(outfile, "\n");
232
/* Close all dpd data files here */
233
for(i=CC_MIN; i <= CC_MAX; i++) psio_close(i,1);
242
char *prgid = "CCLAMBDA";
247
/* put copies of L for excited states in LAMPS with irrep and index label */
248
void Lsave_index(int L_irr, int root_L_irr) {
250
dpdbuf4 L2, LIjAb, LIjbA;
251
char L1A_lbl[32], L1B_lbl[32], L2AA_lbl[32], L2BB_lbl[32], L2AB_lbl[32];
254
sprintf(L1A_lbl, "LIA %d %d", L_irr, root_L_irr);
255
sprintf(L1B_lbl, "Lia %d %d", L_irr, root_L_irr);
256
sprintf(L2AA_lbl, "LIJAB %d %d", L_irr, root_L_irr);
257
sprintf(L2BB_lbl, "Lijab %d %d", L_irr, root_L_irr);
258
sprintf(L2AB_lbl, "LIjAb %d %d", L_irr, root_L_irr);
260
if(params.ref == 0 || params.ref == 1) { /** RHF/ROHF **/
261
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "LIA");
262
dpd_file2_copy(&L1, CC_OEI, L1A_lbl);
263
dpd_file2_close(&L1);
264
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "Lia");
265
dpd_file2_copy(&L1, CC_OEI, L1B_lbl);
266
dpd_file2_close(&L1);
267
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
268
dpd_buf4_copy(&L2, CC_LAMPS, L2AA_lbl);
270
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "Lijab");
271
dpd_buf4_copy(&L2, CC_LAMPS, L2BB_lbl);
273
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
274
dpd_buf4_copy(&L2, CC_LAMPS, L2AB_lbl);
277
else if(params.ref == 2) { /** UHF **/
278
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "LIA");
279
dpd_file2_copy(&L1, CC_OEI, L1A_lbl);
280
dpd_file2_close(&L1);
281
dpd_file2_init(&L1, CC_OEI, L_irr, 2, 3, "Lia");
282
dpd_file2_copy(&L1, CC_OEI, L1B_lbl);
283
dpd_file2_close(&L1);
284
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
285
dpd_buf4_copy(&L2, CC_LAMPS, L2AA_lbl);
287
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "Lijab");
288
dpd_buf4_copy(&L2, CC_LAMPS, L2BB_lbl);
290
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
291
dpd_buf4_copy(&L2, CC_LAMPS, L2AB_lbl);
295
if (params.ref == 0) { /** RHF for those codes that can use them **/
296
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, L2AB_lbl);
297
dpd_buf4_sort(&LIjAb, CC_TMP, pqsr, 0, 5, "LIjbA");
298
sprintf(lbl, "2LIjAb - LIjbA %d %d", L_irr, root_L_irr);
299
dpd_buf4_copy(&LIjAb, CC_LAMPS, lbl);
300
dpd_buf4_close(&LIjAb);
302
sprintf(lbl, "2LIjAb - LIjbA %d %d", L_irr, root_L_irr);
303
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, lbl);
304
dpd_buf4_scm(&LIjAb, 2.0);
305
dpd_buf4_init(&LIjbA, CC_TMP, L_irr, 0, 5, 0, 5, 0, "LIjbA");
306
dpd_buf4_axpy(&LIjbA, &LIjAb, -1.0);
307
dpd_buf4_close(&LIjbA);
308
dpd_buf4_close(&LIjAb);
311
/* also put copy of L1 in LAMPS - RAK thinks this is better organization */
312
if(params.ref == 0 || params.ref == 1) { /** RHF/ROHF **/
313
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "LIA");
314
dpd_file2_copy(&L1, CC_LAMPS, L1A_lbl);
315
dpd_file2_close(&L1);
316
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "Lia");
317
dpd_file2_copy(&L1, CC_LAMPS, L1B_lbl);
318
dpd_file2_close(&L1);
320
else if(params.ref == 2) { /** UHF **/
321
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "LIA");
322
dpd_file2_copy(&L1, CC_LAMPS, L1A_lbl);
323
dpd_file2_close(&L1);
324
dpd_file2_init(&L1, CC_OEI, L_irr, 2, 3, "Lia");
325
dpd_file2_copy(&L1, CC_LAMPS, L1B_lbl);
326
dpd_file2_close(&L1);
331
/* just used for ground state */
332
void L_to_LAMPS(int L_irr) {
334
dpdbuf4 L2, LIjAb, LIjbA;
336
if(params.ref == 0 || params.ref == 1) { /** RHF/ROHF **/
337
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "LIA");
338
dpd_file2_copy(&L1, CC_LAMPS, "LIA");
339
dpd_file2_close(&L1);
340
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "Lia");
341
dpd_file2_copy(&L1, CC_LAMPS, "Lia");
342
dpd_file2_close(&L1);
343
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
344
dpd_buf4_copy(&L2, CC_LAMPS, "LIJAB");
346
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "Lijab");
347
dpd_buf4_copy(&L2, CC_LAMPS, "Lijab");
349
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
350
dpd_buf4_copy(&L2, CC_LAMPS, "LIjAb");
353
else if(params.ref == 2) { /** UHF **/
354
dpd_file2_init(&L1, CC_OEI, L_irr, 0, 1, "LIA");
355
dpd_file2_copy(&L1, CC_OEI, "LIA");
356
dpd_file2_close(&L1);
357
dpd_file2_init(&L1, CC_OEI, L_irr, 2, 3, "Lia");
358
dpd_file2_copy(&L1, CC_OEI, "Lia");
359
dpd_file2_close(&L1);
360
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
361
dpd_buf4_copy(&L2, CC_LAMPS, "LIJAB");
363
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "Lijab");
364
dpd_buf4_copy(&L2, CC_LAMPS, "Lijab");
366
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
367
dpd_buf4_copy(&L2, CC_LAMPS, "LIjAb");
370
if (params.ref == 0) { /** RHF for those codes that can use them **/
371
dpd_buf4_init(&L2, CC_LAMPS, 0, 0, 5, 0, 5, 0, "LIjAb");
372
dpd_buf4_scmcopy(&L2, CC_LAMPS, "2 LIjAb - LIjBa", 2);
373
dpd_buf4_sort_axpy(&L2, CC_LAMPS, pqsr, 0, 5, "2 LIjAb - LIjBa", -1);
377
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, "LIjAb");
378
dpd_buf4_sort(&LIjAb, CC_TMP, pqsr, 0, 5, "LIjbA");
379
dpd_buf4_copy(&LIjAb, CC_LAMPS, "2LIjAb - LIjbA");
380
dpd_buf4_close(&LIjAb);
381
dpd_buf4_init(&LIjAb, CC_LAMPS, L_irr, 0, 5, 0, 5, 0, "2LIjAb - LIjbA");
382
dpd_buf4_scm(&LIjAb, 2.0);
383
dpd_buf4_init(&LIjbA, CC_TMP, L_irr, 0, 5, 0, 5, 0, "LIjbA");
384
dpd_buf4_axpy(&LIjbA, &LIjAb, -1.0);
385
dpd_buf4_close(&LIjbA);
386
dpd_buf4_close(&LIjAb);
392
void L_zero(int L_irr) {
394
dpdbuf4 LIJAB, Lijab, LIjAb;
396
if(params.ref == 0 || params.ref == 1) { /** RHF/ROHF **/
397
dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "New LIA");
398
dpd_file2_init(&Lia, CC_OEI, L_irr, 0, 1, "New Lia");
400
else if(params.ref == 2) { /** UHF **/
401
dpd_file2_init(&LIA, CC_OEI, L_irr, 0, 1, "New LIA");
402
dpd_file2_init(&Lia, CC_OEI, L_irr, 2, 3, "New Lia");
404
dpd_file2_scm(&LIA, 0.0);
405
dpd_file2_scm(&Lia, 0.0);
406
dpd_file2_close(&LIA);
407
dpd_file2_close(&Lia);
409
if (params.ref == 0 || params.ref == 1 ) { /** RHF/ROHF **/
410
dpd_buf4_init(&LIJAB, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
411
dpd_buf4_scm(&LIJAB, 0.0);
412
dpd_buf4_close(&LIJAB);
413
dpd_buf4_init(&Lijab, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New Lijab");
414
dpd_buf4_scm(&Lijab, 0.0);
415
dpd_buf4_close(&Lijab);
416
dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
417
dpd_buf4_scm(&LIjAb, 0.0);
418
dpd_buf4_close(&LIjAb);
421
dpd_buf4_init(&LIJAB, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
422
dpd_buf4_scm(&LIJAB, 0.0);
423
dpd_buf4_close(&LIJAB);
424
dpd_buf4_init(&Lijab, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "New Lijab");
425
dpd_buf4_scm(&Lijab, 0.0);
426
dpd_buf4_close(&Lijab);
427
dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
428
dpd_buf4_scm(&LIjAb, 0.0);
429
dpd_buf4_close(&LIjAb);
added
removed
src/bin/cclambda/cclambda.c
