| [0b990d] | 1 | #
 | 
|---|
 | 2 | eval 'exec perl $0 $*'
 | 
|---|
 | 3 |     if 0;
 | 
|---|
 | 4 | 
 | 
|---|
 | 5 | require QCParse;
 | 
|---|
 | 6 | require Molecule;
 | 
|---|
 | 7 | 
 | 
|---|
 | 8 | ##########################################################################
 | 
|---|
 | 9 | 
 | 
|---|
 | 10 | package QCResult;
 | 
|---|
 | 11 | 
 | 
|---|
 | 12 | # this seems to not work as expected
 | 
|---|
 | 13 | $fltrx = "([-+]?(?:\\d+\\.\\d*|\\d+|\\.\\d+)(?:[eEdD][+-]?\\d+)?)";
 | 
|---|
 | 14 | 
 | 
|---|
 | 15 | $have_nodenum = 0;
 | 
|---|
 | 16 | 
 | 
|---|
 | 17 | sub test {
 | 
|---|
 | 18 |     my $i;
 | 
|---|
 | 19 |     @nums = ("-1.0", "-1", "-.1", "-1.",
 | 
|---|
 | 20 |              "-1.0E-04", "-1E-04", "-.1E-04", "-1.E-04");
 | 
|---|
 | 21 |     foreach $i ( @nums ) {
 | 
|---|
 | 22 |         $i =~ /$fltrx/;
 | 
|---|
 | 23 |         my $flt = $1;
 | 
|---|
 | 24 |         printf "%10s %10s %12.8f\n", $i, $flt, $flt;
 | 
|---|
 | 25 |     }
 | 
|---|
 | 26 | }
 | 
|---|
 | 27 | 
 | 
|---|
 | 28 | sub new {
 | 
|---|
 | 29 |     my $this = shift;
 | 
|---|
 | 30 |     my $class = ref($this) || $this;
 | 
|---|
 | 31 |     my $self = {};
 | 
|---|
 | 32 |     bless $self, $class;
 | 
|---|
 | 33 |     $self->initialize(@_);
 | 
|---|
 | 34 |     return $self;
 | 
|---|
 | 35 | }
 | 
|---|
 | 36 | 
 | 
|---|
 | 37 | sub initialize {
 | 
|---|
 | 38 |     my $self = shift;
 | 
|---|
 | 39 |     my $infile = shift;
 | 
|---|
 | 40 |     my $outfile = shift;
 | 
|---|
 | 41 |     my $parse = new QCParse();
 | 
|---|
 | 42 |     $parse->parse_file($infile);
 | 
|---|
 | 43 |     $self->{"qcinput"} = QCInput::new QCInput($parse);
 | 
|---|
 | 44 | 
 | 
|---|
 | 45 |     $self->{"exists"} = 0;
 | 
|---|
 | 46 |     $self->{"ok"} = 0;
 | 
|---|
 | 47 |     if (-e $outfile) {
 | 
|---|
 | 48 |         $self->{"exists"} = 1;
 | 
|---|
 | 49 |         open(OUTFILE,"<$outfile");
 | 
|---|
 | 50 |         my $first = 1;
 | 
|---|
 | 51 |         while (<OUTFILE>) {
 | 
|---|
 | 52 |             if ($first && /^ *[0-9]+:/) {
 | 
|---|
 | 53 |                 $have_nodenum = 1;
 | 
|---|
 | 54 |                 $first = 0;
 | 
|---|
 | 55 |             }
 | 
|---|
 | 56 |             s/^ *[0-9]+:// if ($have_nodenum);
 | 
|---|
 | 57 |             if (/^\s*MPQC:/) {
 | 
|---|
 | 58 |                 $self->parse_mpqc(\*OUTFILE);
 | 
|---|
 | 59 |                 last;
 | 
|---|
 | 60 |             }
 | 
|---|
 | 61 |             elsif (/^\s*Entering Gaussian System/) {
 | 
|---|
 | 62 |                 $self->parse_g94(\*OUTFILE);
 | 
|---|
 | 63 |                 last;
 | 
|---|
 | 64 |             }
 | 
|---|
 | 65 |         }
 | 
|---|
 | 66 |         close(OUTFILE);
 | 
|---|
 | 67 |     }
 | 
|---|
 | 68 | }
 | 
|---|
 | 69 | 
 | 
|---|
 | 70 | sub parse_g94 {
 | 
|---|
 | 71 |     my $self = shift;
 | 
|---|
 | 72 |     my $out = shift;
 | 
|---|
 | 73 |     my $scfenergy = "";
 | 
|---|
 | 74 |     my $mp2energy = "";
 | 
|---|
 | 75 |     my $ccsd_tenergy = "";
 | 
|---|
 | 76 |     my $t1norm = "";
 | 
|---|
 | 77 |     my $optconverged = 0;
 | 
|---|
 | 78 |     my $molecule;
 | 
|---|
 | 79 |     my $havefreq = 0;
 | 
|---|
 | 80 |     my $freq = [];
 | 
|---|
 | 81 |     my $ifreq = 0;
 | 
|---|
 | 82 |     my $b3pw91energy = "";
 | 
|---|
 | 83 |     while (<$out>) {
 | 
|---|
 | 84 |         s/^ *[0-9]+:// if ($have_nodenum);
 | 
|---|
 | 85 |         if (/^\s*SCF Done:  E\(RHF\) =\s*$fltrx\s/) {
 | 
|---|
 | 86 |             $scfenergy = $1;
 | 
|---|
 | 87 |         }
 | 
|---|
 | 88 |         elsif (/^\s*E2\s*=\s*$fltrx\s*EUMP2\s*=\s*$fltrx/) {
 | 
|---|
 | 89 |             $mp2energy = $2;
 | 
|---|
 | 90 |         }
 | 
|---|
 | 91 |         elsif (/^\s*CCSD\(T\)\s*=\s*$fltrx/) {
 | 
|---|
 | 92 |             $ccsd_tenergy = $1;
 | 
|---|
 | 93 |             $ccsd_tenergy =~ s/[DdE]/e/;
 | 
|---|
 | 94 |         }
 | 
|---|
 | 95 |         elsif (/E\(RB\+HF-PW91\)\s*=\s*$fltrx/) {
 | 
|---|
 | 96 |             $b3pw91energy = $1;
 | 
|---|
 | 97 |             $b3pw91energy =~ s/[DdE]/e/;
 | 
|---|
 | 98 |         }
 | 
|---|
 | 99 |         elsif (/^\s*T1 Diagnostic\s*=\s*$fltrx/) {
 | 
|---|
 | 100 |             $t1norm = $1;
 | 
|---|
 | 101 |         }
 | 
|---|
 | 102 |         elsif (/^\s*-- Stationary point found./
 | 
|---|
 | 103 |                || /CONVERGENCE CRITERIA APPARENTLY SATISFIED/) {
 | 
|---|
 | 104 |             $optconverged = 1;
 | 
|---|
 | 105 |         }
 | 
|---|
 | 106 |         elsif ($optconverged
 | 
|---|
 | 107 |                && (/^\s*Input orientation:/ || /^\s*Standard orientation:/)) {
 | 
|---|
 | 108 |             <$out>; <$out>; <$out>; <$out>;
 | 
|---|
 | 109 |             my $molstr = "";
 | 
|---|
 | 110 |             while (<$out>) {
 | 
|---|
 | 111 |                 if (! /^\s+\d+\s+(\d+)\s+$fltrx\s+$fltrx\s+$fltrx\s+/) {
 | 
|---|
 | 112 |                     last;
 | 
|---|
 | 113 |                 }
 | 
|---|
 | 114 |                 $molstr = "${molstr} $1 $2 $3 $4\n";
 | 
|---|
 | 115 |             }
 | 
|---|
 | 116 |             $molecule = new Molecule($molstr);
 | 
|---|
 | 117 |         }
 | 
|---|
 | 118 |         elsif (/^\s*Frequencies --\s+$fltrx\s+$fltrx\s+$fltrx/) {
 | 
|---|
 | 119 |             $freq->[$ifreq] = $1; $ifreq++;
 | 
|---|
 | 120 |             $freq->[$ifreq] = $2; $ifreq++;
 | 
|---|
 | 121 |             $freq->[$ifreq] = $3; $ifreq++;
 | 
|---|
 | 122 |             $havefreq = 1;
 | 
|---|
 | 123 |         }
 | 
|---|
 | 124 |         elsif (/^\s*Frequencies --\s+$fltrx\s+$fltrx/) {
 | 
|---|
 | 125 |             $freq->[$ifreq] = $1; $ifreq++;
 | 
|---|
 | 126 |             $freq->[$ifreq] = $2; $ifreq++;
 | 
|---|
 | 127 |             $havefreq = 1;
 | 
|---|
 | 128 |         }
 | 
|---|
 | 129 |         elsif (/^\s*Frequencies --\s+$fltrx\s/) {
 | 
|---|
 | 130 |             $freq->[$ifreq] = $1; $ifreq++;
 | 
|---|
 | 131 |             $havefreq = 1;
 | 
|---|
 | 132 |         }
 | 
|---|
 | 133 |     }
 | 
|---|
 | 134 |     $self->{"scfenergy"} = $scfenergy;
 | 
|---|
 | 135 |     $self->{"mp2energy"} = $mp2energy;
 | 
|---|
 | 136 |     $self->{"b3pw91energy"} = $b3pw91energy;
 | 
|---|
 | 137 |     $self->{"optconverged"} = $optconverged;
 | 
|---|
 | 138 |     if ($optconverged) {
 | 
|---|
 | 139 |         $self->{"optmolecule"} = $molecule;
 | 
|---|
 | 140 |     }
 | 
|---|
 | 141 |     $self->{"have_frequencies"} = $havefreq;
 | 
|---|
 | 142 |     if ($havefreq) {
 | 
|---|
 | 143 |         my @tmp = sort(@$freq);
 | 
|---|
 | 144 |         $freq = \@tmp;
 | 
|---|
 | 145 |         $self->{"freq"} = $freq;
 | 
|---|
 | 146 |     }
 | 
|---|
 | 147 |     my $qcinput = $self->{"qcinput"};
 | 
|---|
 | 148 |     my $method = $qcinput->method();
 | 
|---|
 | 149 |     if ($method eq "MP2") {
 | 
|---|
 | 150 |         $self->{"energy"} = $mp2energy;
 | 
|---|
 | 151 |     }
 | 
|---|
 | 152 |     elsif ($method eq "CCSD(T)") {
 | 
|---|
 | 153 |         $self->{"energy"} = $ccsd_tenergy;
 | 
|---|
 | 154 |     }
 | 
|---|
 | 155 |     elsif ($method eq "B3-PW91") {
 | 
|---|
 | 156 |         $self->{"energy"} = $b3pw91energy;
 | 
|---|
 | 157 |     }
 | 
|---|
 | 158 |     elsif ($method eq "SCF"
 | 
|---|
 | 159 |            || $method eq "ROSCF") {
 | 
|---|
 | 160 |         $self->{"energy"} = $scfenergy;
 | 
|---|
 | 161 |     }
 | 
|---|
 | 162 | 
 | 
|---|
 | 163 |     $self->{"t1norm"} = $t1norm;
 | 
|---|
 | 164 | 
 | 
|---|
 | 165 |     $self->{"ok"} = 0;
 | 
|---|
 | 166 |     if ($self->{"energy"} ne "") {
 | 
|---|
 | 167 |         if ($qcinput->optimize()) {
 | 
|---|
 | 168 |             if ($self->{"optconverged"}) {
 | 
|---|
 | 169 |                 $self->{"ok"} = 1
 | 
|---|
 | 170 |                 }
 | 
|---|
 | 171 |             else {
 | 
|---|
 | 172 |                 #printf "not ok because not converged\n";
 | 
|---|
 | 173 |             }
 | 
|---|
 | 174 |         }
 | 
|---|
 | 175 |         else {
 | 
|---|
 | 176 |             $self->{"ok"} = 1;
 | 
|---|
 | 177 |         }
 | 
|---|
 | 178 |         if ($qcinput->frequencies() && ! $havefreq) {
 | 
|---|
 | 179 |             $self->{"ok"} = 0;
 | 
|---|
 | 180 |             #printf "not ok because no freq\n";
 | 
|---|
 | 181 |         }
 | 
|---|
 | 182 |     }
 | 
|---|
 | 183 |     else {
 | 
|---|
 | 184 |         #printf "not ok because no energy\n";
 | 
|---|
 | 185 |     }
 | 
|---|
 | 186 | }
 | 
|---|
 | 187 | 
 | 
|---|
 | 188 | sub parse_mpqc {
 | 
|---|
 | 189 |     my $self = shift;
 | 
|---|
 | 190 |     my $out = shift;
 | 
|---|
 | 191 |     my $scfenergy = "";
 | 
|---|
 | 192 |     my $opt1energy = "";
 | 
|---|
 | 193 |     my $opt2energy = "";
 | 
|---|
 | 194 |     my $zapt2energy = "";
 | 
|---|
 | 195 |     my $mp2energy = "";
 | 
|---|
 | 196 |     my $optconverged = 0;
 | 
|---|
 | 197 |     my $molecule;
 | 
|---|
 | 198 |     my $havefreq = 0;
 | 
|---|
 | 199 |     my $freq;
 | 
|---|
 | 200 |     my $wante = 1;
 | 
|---|
 | 201 |     my $ifreq = 0;
 | 
|---|
 | 202 |     my $grad;
 | 
|---|
 | 203 |     my $ngrad;
 | 
|---|
 | 204 |     my $state = "none";
 | 
|---|
 | 205 |     my $molecularenergy = "";
 | 
|---|
 | 206 |     my $s2norm = "";
 | 
|---|
 | 207 |     my $to_angstrom = 0.52917706;
 | 
|---|
 | 208 |     my $geometry_conversion = $to_angstrom;
 | 
|---|
 | 209 |     my $error = "";
 | 
|---|
 | 210 |     my $ccsd_energy = "";
 | 
|---|
 | 211 |     my $ccsd_t_energy = "";
 | 
|---|
 | 212 |     my $t1norm = "";
 | 
|---|
 | 213 |     my $t12norm = "";
 | 
|---|
 | 214 |     my $psioutput = "";
 | 
|---|
 | 215 |     while (<$out>) {
 | 
|---|
 | 216 |         s/^ *[0-9]+:// if ($have_nodenum);
 | 
|---|
 | 217 |         if ($state eq "read grad" && $wante) {
 | 
|---|
 | 218 |             if (/^\s*([0-9]+\s+[A-Za-z]+\s+)?([\-\.0-9]+)\s+([\-\.0-9]+)\s+([\-\.0-9]+)/) {
 | 
|---|
 | 219 |                 $grad->[$ngrad + 0] = $2;
 | 
|---|
 | 220 |                 $grad->[$ngrad + 1] = $3;
 | 
|---|
 | 221 |                 $grad->[$ngrad + 2] = $4;
 | 
|---|
 | 222 |                 $ngrad = $ngrad + 3;
 | 
|---|
 | 223 |             }
 | 
|---|
 | 224 |             elsif (/^\s*([0-9]+\s+)?([\-\.0-9]+)/) {
 | 
|---|
 | 225 |                 $grad->[$ngrad] = $2;
 | 
|---|
 | 226 |                 $ngrad = $ngrad + 1;
 | 
|---|
 | 227 |             }
 | 
|---|
 | 228 |             else {
 | 
|---|
 | 229 |                 $self->{"grad"} = $grad;
 | 
|---|
 | 230 |                 $state = "none";
 | 
|---|
 | 231 |             }
 | 
|---|
 | 232 |         }
 | 
|---|
 | 233 | 
 | 
|---|
 | 234 |         if ($wante && /total scf energy =\s+$fltrx/) {
 | 
|---|
 | 235 |             $scfenergy = $1;
 | 
|---|
 | 236 |         }
 | 
|---|
 | 237 |         elsif ($wante && /OPT1 energy .*:\s+$fltrx/) {
 | 
|---|
 | 238 |             $opt1energy = $1;
 | 
|---|
 | 239 |         }
 | 
|---|
 | 240 |         elsif ($wante && /OPT2 energy .*:\s+$fltrx/) {
 | 
|---|
 | 241 |             $opt2energy = $1;
 | 
|---|
 | 242 |         }
 | 
|---|
 | 243 |         elsif ($wante && /ZAPT2 energy .*:\s+$fltrx/) {
 | 
|---|
 | 244 |             $zapt2energy = $1;
 | 
|---|
 | 245 |         }
 | 
|---|
 | 246 |         elsif ($wante && /MP2 energy .*:\s+$fltrx/) {
 | 
|---|
 | 247 |             $mp2energy = $1;
 | 
|---|
 | 248 |         }
 | 
|---|
 | 249 |         elsif (/CSCF: An SCF program written in C/) {
 | 
|---|
 | 250 |             $psioutput = 1;
 | 
|---|
 | 251 |         }
 | 
|---|
 | 252 |         elsif ($psioutput && $wante
 | 
|---|
 | 253 |                && /total energy       =\s*$fltrx/) {
 | 
|---|
 | 254 |             $scfenergy = $1;
 | 
|---|
 | 255 |         }
 | 
|---|
 | 256 |         elsif ($psioutput && $wante
 | 
|---|
 | 257 |                && /\s*ITER\s+CORRELATION ENERGY\s+T1 2-NORM\s+T1 DIAG/) {
 | 
|---|
 | 258 |             # this is a PSI CC output embedded in an MPQC output
 | 
|---|
 | 259 |             # grab iteration 0
 | 
|---|
 | 260 |             my $ecorr;
 | 
|---|
 | 261 |             <$out>;
 | 
|---|
 | 262 |             while (<$out>) {
 | 
|---|
 | 263 |                 if (/^\s*\d+\s+$fltrx\s+$fltrx\s+$fltrx/) {
 | 
|---|
 | 264 |                     $ecorr = $1;
 | 
|---|
 | 265 |                     $t12norm = $2;
 | 
|---|
 | 266 |                     $t1norm = $3;
 | 
|---|
 | 267 |                 }
 | 
|---|
 | 268 |                 else {
 | 
|---|
 | 269 |                     last;
 | 
|---|
 | 270 |                 }
 | 
|---|
 | 271 |             }
 | 
|---|
 | 272 |             $ccsd_energy = $ecorr + $scfenergy;
 | 
|---|
 | 273 |         }
 | 
|---|
 | 274 | # IMBN: Added the following elsif to handle new cc output 
 | 
|---|
 | 275 |         elsif ($psioutput && $wante
 | 
|---|
 | 276 |                && /\s*ITER\s+CORRELATION ENERGY\s+T1 DIAG\s+D1\(CCSD\)/) {
 | 
|---|
 | 277 |             # this is a PSI CC output embedded in an MPQC output
 | 
|---|
 | 278 |             # grab iteration 0
 | 
|---|
 | 279 |             my $ecorr;
 | 
|---|
 | 280 |             <$out>;
 | 
|---|
 | 281 |             while (<$out>) {
 | 
|---|
 | 282 |                 if (/^\s*\d+\s+$fltrx\s+$fltrx\s+$fltrx/) {
 | 
|---|
 | 283 |                     $ecorr = $1;
 | 
|---|
 | 284 |                     $t1norm = $2;
 | 
|---|
 | 285 |                     $t12norm = $3;
 | 
|---|
 | 286 |                 }
 | 
|---|
 | 287 |                 else {
 | 
|---|
 | 288 |                     last;
 | 
|---|
 | 289 |                 }
 | 
|---|
 | 290 |             }
 | 
|---|
 | 291 |             $ccsd_energy = $ecorr + $scfenergy;
 | 
|---|
 | 292 |         }
 | 
|---|
 | 293 |         elsif ($wante && /Value of the MolecularEnergy:\s+$fltrx/) {
 | 
|---|
 | 294 |             $molecularenergy = $1;
 | 
|---|
 | 295 |         }
 | 
|---|
 | 296 |         elsif ($wante && /S2\(ov\) matrix 1-norm\s*=\s*$fltrx/) {
 | 
|---|
 | 297 |             $self->{"s2matrix1norm"} = $1;
 | 
|---|
 | 298 |         }
 | 
|---|
 | 299 |         elsif ($wante && /There are degenerate orbitals within an irrep/) {
 | 
|---|
 | 300 |             $self->{"degenerate"} = 1;
 | 
|---|
 | 301 |         }
 | 
|---|
 | 302 |         elsif ($wante && /D1\(MP2\)\s*=\s*$fltrx/) {
 | 
|---|
 | 303 |             $self->{"d1mp2"} = $1;
 | 
|---|
 | 304 |         }
 | 
|---|
 | 305 |         elsif ($wante && /D2\(MP1\)\s*=\s*$fltrx/) {
 | 
|---|
 | 306 |             $self->{"d2mp1"} = $1;
 | 
|---|
 | 307 |         }
 | 
|---|
 | 308 |         elsif ($wante && /S2\(ov\) matrix inf-norm\s*=\s*$fltrx/) {
 | 
|---|
 | 309 |             $self->{"s2matrixinfnorm"} = $1;
 | 
|---|
 | 310 |         }
 | 
|---|
 | 311 |         elsif ($wante && /S2 norm\s*=\s*$fltrx/) {
 | 
|---|
 | 312 |             $s2norm = $1;
 | 
|---|
 | 313 |         }
 | 
|---|
 | 314 |         elsif ($wante && /Largest S2 values.*:/) {
 | 
|---|
 | 315 |             my $s2large_coef = [];
 | 
|---|
 | 316 |             my $s2large_i = [];
 | 
|---|
 | 317 |             my $s2large_a = [];
 | 
|---|
 | 318 |             my $is2 = 0;
 | 
|---|
 | 319 |             while (<$out>) {
 | 
|---|
 | 320 |                 if (/^\s*([0-9]+)\s+$fltrx\s+([0-9]+)\s+([A-Za-z0-9\'\"]+)\s+->\s+([0-9]+)\s+([A-Za-z0-9\'\"]+)\s*$/) {
 | 
|---|
 | 321 |                     $s2large_coef->[$is2] = $2;
 | 
|---|
 | 322 |                     $s2large_i->[$is2] = "$3$4";
 | 
|---|
 | 323 |                     $s2large_a->[$is2] = "$5$6";
 | 
|---|
 | 324 |                     $is2 = $is2 + 1;
 | 
|---|
 | 325 |                 }
 | 
|---|
 | 326 |                 else {
 | 
|---|
 | 327 |                     last;
 | 
|---|
 | 328 |                 }
 | 
|---|
 | 329 |             }
 | 
|---|
 | 330 |             $self->{"s2large_coef"} = $s2large_coef;
 | 
|---|
 | 331 |             $self->{"s2large_i"} = $s2large_i;
 | 
|---|
 | 332 |             $self->{"s2large_a"} = $s2large_a;
 | 
|---|
 | 333 |         }
 | 
|---|
 | 334 |         elsif ($wante && /Largest first order coefficients.*:/) {
 | 
|---|
 | 335 |             my $d1large_coef = [];
 | 
|---|
 | 336 |             my $d1large_i = [];
 | 
|---|
 | 337 |             my $d1large_j = [];
 | 
|---|
 | 338 |             my $d1large_a = [];
 | 
|---|
 | 339 |             my $d1large_b = [];
 | 
|---|
 | 340 |             my $d1large_spin = [];
 | 
|---|
 | 341 |             my $id1 = 0;
 | 
|---|
 | 342 |             while (<$out>) {
 | 
|---|
 | 343 |                 if (/^\s*([0-9]+)\s+$fltrx\s+([0-9]+)\s+([A-Za-z0-9\'\"]+)\s+([0-9]+)\s+([A-Za-z0-9\'\"]+)\s+->\s+([0-9]+)\s+([A-Za-z0-9\'\"]+)\s+([0-9]+)\s+([A-Za-z0-9\'\"]+)\s+\((....)\)\s*$/) {
 | 
|---|
 | 344 |                     $d1large_coef->[$id1] = $2;
 | 
|---|
 | 345 |                     $d1large_i->[$id1] = "$3$4";
 | 
|---|
 | 346 |                     $d1large_j->[$id1] = "$5$6";
 | 
|---|
 | 347 |                     $d1large_a->[$id1] = "$7$8";
 | 
|---|
 | 348 |                     $d1large_b->[$id1] = "$9$10";
 | 
|---|
 | 349 |                     $d1large_spin->[$id1] = $11;
 | 
|---|
 | 350 |                     $id1 = $id1 + 1;
 | 
|---|
 | 351 |                 }
 | 
|---|
 | 352 |                 else {
 | 
|---|
 | 353 |                     last;
 | 
|---|
 | 354 |                 }
 | 
|---|
 | 355 |             }
 | 
|---|
 | 356 |             $self->{"d1large_coef"} = $d1large_coef;
 | 
|---|
 | 357 |             $self->{"d1large_i"} = $d1large_i;
 | 
|---|
 | 358 |             $self->{"d1large_j"} = $d1large_j;
 | 
|---|
 | 359 |             $self->{"d1large_a"} = $d1large_a;
 | 
|---|
 | 360 |             $self->{"d1large_b"} = $d1large_b;
 | 
|---|
 | 361 |             $self->{"d1large_spin"} = $d1large_spin;
 | 
|---|
 | 362 |         }
 | 
|---|
 | 363 |         elsif ($wante && /^\s*Natural\s+Population\s+Analysis:\s*$/) {
 | 
|---|
 | 364 |             my $npacharge = [];
 | 
|---|
 | 365 |             my $npashellpop = [];
 | 
|---|
 | 366 |             <$out>;
 | 
|---|
 | 367 |             my $iatom = 0;
 | 
|---|
 | 368 |             while (<$out>) {
 | 
|---|
 | 369 |                 if (/^\s*\d+\s+[A-Za-z]+\s+$fltrx\s+$fltrx\s*$/) {
 | 
|---|
 | 370 |                     $npacharge->[$iatom] = $1;
 | 
|---|
 | 371 |                     $npashellpop->[$iatom] = [ $2 ];
 | 
|---|
 | 372 |                 }
 | 
|---|
 | 373 |                 elsif (/^\s*\d+\s+[A-Za-z]+\s+$fltrx\s+$fltrx\s+$fltrx\s*$/) {
 | 
|---|
 | 374 |                     $npacharge->[$iatom] = $1;
 | 
|---|
 | 375 |                     $npashellpop->[$iatom] = [ $2, $3 ];
 | 
|---|
 | 376 |                 }
 | 
|---|
 | 377 |                 elsif (/^\s*\d+\s+[A-Za-z]+\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s*$/) {
 | 
|---|
 | 378 |                     $npacharge->[$iatom] = $1;
 | 
|---|
 | 379 |                     $npashellpop->[$iatom] = [ $2, $3, $4 ];
 | 
|---|
 | 380 |                 }
 | 
|---|
 | 381 |                 elsif (/^\s*\d+\s+[A-Za-z]+\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s*$/) {
 | 
|---|
 | 382 |                     $npacharge->[$iatom] = $1;
 | 
|---|
 | 383 |                     $npashellpop->[$iatom] = [ $2, $3, $4, $5 ];
 | 
|---|
 | 384 |                 }
 | 
|---|
 | 385 |                 elsif (/^\s*\d+\s+[A-Za-z]+\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s*$/) {
 | 
|---|
 | 386 |                     $npacharge->[$iatom] = $1;
 | 
|---|
 | 387 |                     $npashellpop->[$iatom] = [ $2, $3, $4, $5, $6 ];
 | 
|---|
 | 388 |                 }
 | 
|---|
 | 389 |                 elsif (/^\s*\d+\s+[A-Za-z]+\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s+$fltrx\s*$/) {
 | 
|---|
 | 390 |                     $npacharge->[$iatom] = $1;
 | 
|---|
 | 391 |                     $npashellpop->[$iatom] = [ $2, $3, $4, $5, $6, $7 ];
 | 
|---|
 | 392 |                 }
 | 
|---|
 | 393 |                 elsif (/^\s*$/) {
 | 
|---|
 | 394 |                     last;
 | 
|---|
 | 395 |                 }
 | 
|---|
 | 396 |                 else {
 | 
|---|
 | 397 |                     die "AM too high to read NPA shell populations (line=$_)";
 | 
|---|
 | 398 |                 }
 | 
|---|
 | 399 |                 $iatom = $iatom + 1;
 | 
|---|
 | 400 |             }
 | 
|---|
 | 401 |             $self->{"npacharge"} = $npacharge;
 | 
|---|
 | 402 |             $self->{"npashellpop"} = $npashellpop;
 | 
|---|
 | 403 |         }
 | 
|---|
 | 404 |         elsif (/The optimization has converged/) {
 | 
|---|
 | 405 |             $optconverged = 1;
 | 
|---|
 | 406 |         }
 | 
|---|
 | 407 |         elsif (/^\s*Beginning displacement/) {
 | 
|---|
 | 408 |             # don't grab energies for displaced geometries
 | 
|---|
 | 409 |             $wante = 0;
 | 
|---|
 | 410 |         }
 | 
|---|
 | 411 |         elsif ($optconverged
 | 
|---|
 | 412 |                && /^\s+n\s+atom\s+label\s+x\s+y\s+z\s+mass\s*$/) {
 | 
|---|
 | 413 |             # old style geometry
 | 
|---|
 | 414 |             my $molstr = "";
 | 
|---|
 | 415 |             while (<$out>) {
 | 
|---|
 | 416 |                 s/^ *[0-9]+:// if ($have_nodenum);
 | 
|---|
 | 417 |                 if (! /^\s+\d+\s+(\w+)\s+$fltrx\s+$fltrx\s+$fltrx\s+/
 | 
|---|
 | 418 |                  && ! /^\s+\d+\s+(\w+)\s+\S+\s+$fltrx\s+$fltrx\s+$fltrx\s+/) {
 | 
|---|
 | 419 |                     last;
 | 
|---|
 | 420 |                 }
 | 
|---|
 | 421 |                 $molstr = sprintf "%s %s %16.14f %16.14f %16.14f\n",
 | 
|---|
 | 422 |                     ${molstr}, $1, $2 * $to_angstrom,
 | 
|---|
 | 423 |                     $3 * $to_angstrom, $4 * $to_angstrom;
 | 
|---|
 | 424 |             }
 | 
|---|
 | 425 |             $molecule = new Molecule($molstr);
 | 
|---|
 | 426 |         }
 | 
|---|
 | 427 |         elsif (/^\s*molecule<Molecule>:/) {
 | 
|---|
 | 428 |             $geometry_conversion = $to_angstrom;
 | 
|---|
 | 429 |         }
 | 
|---|
 | 430 |         elsif (/^\s*unit = \"angstrom\"\s*/) { # "
 | 
|---|
 | 431 |             $geometry_conversion = 1.0;
 | 
|---|
 | 432 |         }
 | 
|---|
 | 433 |         elsif ($optconverged
 | 
|---|
 | 434 |                && /n\s+atoms\s+(atom_labels\s+)?geometry/) {
 | 
|---|
 | 435 |             # new style geometry
 | 
|---|
 | 436 |             my $molstr = "";
 | 
|---|
 | 437 |             while (<$out>) {
 | 
|---|
 | 438 |                 s/^ *[0-9]+:// if ($have_nodenum);
 | 
|---|
 | 439 |                 if (! /^\s+\d+\s+(\w+)\s+\[\s*$fltrx\s+$fltrx\s+$fltrx\s*\]/
 | 
|---|
 | 440 |                  && ! /^\s+\d+\s+(\w+)\s+\"[^\"]*\"+\s+\[\s*$fltrx\s+$fltrx\s+$fltrx\s*\]/) { # " (unconfuse emacs)
 | 
|---|
 | 441 |                     last;
 | 
|---|
 | 442 |                 }
 | 
|---|
 | 443 |                 $molstr = sprintf "%s %s %16.14f %16.14f %16.14f\n",
 | 
|---|
 | 444 |                     ${molstr}, $1, $2 * $geometry_conversion,
 | 
|---|
 | 445 |                     $3 * $geometry_conversion, $4 * $geometry_conversion;
 | 
|---|
 | 446 |             }
 | 
|---|
 | 447 |             $molecule = new Molecule($molstr);
 | 
|---|
 | 448 |         }
 | 
|---|
 | 449 |         elsif (/^\s+Total (MP2 )?[Gg]radient/) {
 | 
|---|
 | 450 |             $state = "read grad";
 | 
|---|
 | 451 |             $grad = [];
 | 
|---|
 | 452 |             $ngrad = 0;
 | 
|---|
 | 453 |         }
 | 
|---|
 | 454 |         elsif (/^\s+Frequencies .*:\s*$/) {
 | 
|---|
 | 455 |             # read the frequencies
 | 
|---|
 | 456 |             while (<$out>) {
 | 
|---|
 | 457 |                 s/^ *[0-9]+:// if ($have_nodenum);
 | 
|---|
 | 458 |                 if (/^\s+\d+\s+$fltrx\s*$/) {
 | 
|---|
 | 459 |                     $freq->[$ifreq] = $1;
 | 
|---|
 | 460 |                     $ifreq++;
 | 
|---|
 | 461 |                 }
 | 
|---|
 | 462 |                 elsif (/THERMODYNAMIC ANALYSIS/) {
 | 
|---|
 | 463 |                     last;
 | 
|---|
 | 464 |                 }
 | 
|---|
 | 465 |             }
 | 
|---|
 | 466 |             $havefreq = 1;
 | 
|---|
 | 467 |         }
 | 
|---|
 | 468 |     }
 | 
|---|
 | 469 |     $self->{"t1norm"} = $t1norm;
 | 
|---|
 | 470 |     $self->{"t1matrix2norm"} = $t12norm;
 | 
|---|
 | 471 |     $self->{"s2norm"} = $s2norm;
 | 
|---|
 | 472 |     $self->{"ccsdenergy"} = $ccsd_energy;
 | 
|---|
 | 473 |     $self->{"ccsd_tenergy"} = $ccsd_t_energy;
 | 
|---|
 | 474 |     $self->{"scfenergy"} = $scfenergy;
 | 
|---|
 | 475 |     $self->{"opt1energy"} = $opt1energy;
 | 
|---|
 | 476 |     $self->{"opt2energy"} = $opt2energy;
 | 
|---|
 | 477 |     $self->{"zapt2energy"} = $zapt2energy;
 | 
|---|
 | 478 |     $self->{"mp2energy"} = $mp2energy;
 | 
|---|
 | 479 |     $self->{"optconverged"} = $optconverged;
 | 
|---|
 | 480 |     $self->{"molecularenergy"} = $molecularenergy;
 | 
|---|
 | 481 |     if ($optconverged) {
 | 
|---|
 | 482 |         $self->{"optmolecule"} = $molecule;
 | 
|---|
 | 483 |     }
 | 
|---|
 | 484 |     $self->{"have_frequencies"} = $havefreq;
 | 
|---|
 | 485 |     if ($havefreq) {
 | 
|---|
 | 486 |         my @tmp = sort(@$freq);
 | 
|---|
 | 487 |         $freq = \@tmp;
 | 
|---|
 | 488 |         $self->{"freq"} = $freq;
 | 
|---|
 | 489 |     }
 | 
|---|
 | 490 |     my $qcinput = $self->{"qcinput"};
 | 
|---|
 | 491 |     my $method = $qcinput->method();
 | 
|---|
 | 492 |     #printf "qcinput ok = %d\n", $qcinput->ok();
 | 
|---|
 | 493 |     if ($method eq "MP2") {
 | 
|---|
 | 494 |         if ($mp2energy ne "") {
 | 
|---|
 | 495 |             $self->{"energy"} = $mp2energy;
 | 
|---|
 | 496 |         }
 | 
|---|
 | 497 |         elsif ($qcinput->mult() == 1) {
 | 
|---|
 | 498 |             $self->{"energy"} = $zapt2energy;
 | 
|---|
 | 499 |         }
 | 
|---|
 | 500 |     }
 | 
|---|
 | 501 |     elsif ($method eq "OPT1[2]") {
 | 
|---|
 | 502 |         $self->{"energy"} = $opt1energy;
 | 
|---|
 | 503 |     }
 | 
|---|
 | 504 |     elsif ($method eq "OPT2[2]") {
 | 
|---|
 | 505 |         $self->{"energy"} = $opt2energy;
 | 
|---|
 | 506 |     }
 | 
|---|
 | 507 |     elsif ($method eq "CCSD") {
 | 
|---|
 | 508 |         $self->{"energy"} = $ccsd_energy;
 | 
|---|
 | 509 |     }
 | 
|---|
 | 510 |     elsif ($method eq "CCSD(T)") {
 | 
|---|
 | 511 |         $self->{"energy"} = $ccsd_t_energy;
 | 
|---|
 | 512 |     }
 | 
|---|
 | 513 |     elsif ($qcinput->ok()
 | 
|---|
 | 514 |            && ($method eq "SCF" || $method eq "ROSCF")) {
 | 
|---|
 | 515 |         $self->{"energy"} = $scfenergy;
 | 
|---|
 | 516 |     }
 | 
|---|
 | 517 | 
 | 
|---|
 | 518 |     if ($self->{"energy"} eq "") {
 | 
|---|
 | 519 |         $self->{"energy"} = $self->{"molecularenergy"};
 | 
|---|
 | 520 |     }
 | 
|---|
 | 521 | 
 | 
|---|
 | 522 |     $self->{"ok"} = 0;
 | 
|---|
 | 523 |     if ($self->{"energy"} ne "") {
 | 
|---|
 | 524 |         if ($qcinput->optimize()) {
 | 
|---|
 | 525 |             if ($self->{"optconverged"}) {
 | 
|---|
 | 526 |                 $self->{"ok"} = 1;
 | 
|---|
 | 527 |             }
 | 
|---|
 | 528 |             else {
 | 
|---|
 | 529 |                 $error = "$error geom not converged\n";
 | 
|---|
 | 530 |             }
 | 
|---|
 | 531 |         }
 | 
|---|
 | 532 |         else {
 | 
|---|
 | 533 |             $self->{"ok"} = 1;
 | 
|---|
 | 534 |         }
 | 
|---|
 | 535 |         if ($qcinput->frequencies() && ! $havefreq) {
 | 
|---|
 | 536 |             $self->{"ok"} = 0;
 | 
|---|
 | 537 |             $error = "$error no freq\n";
 | 
|---|
 | 538 |         }
 | 
|---|
 | 539 |     }
 | 
|---|
 | 540 |     else {
 | 
|---|
 | 541 |         $error = "$error no energy\n";
 | 
|---|
 | 542 |     }
 | 
|---|
 | 543 | 
 | 
|---|
 | 544 |     if (! $qcinput->ok()) {
 | 
|---|
 | 545 |         $self->{"ok"} = 0;
 | 
|---|
 | 546 |         $error = "$error qcinput error: $qcinput->{'error'}";
 | 
|---|
 | 547 |     }
 | 
|---|
 | 548 | 
 | 
|---|
 | 549 |     $self->{"error"} = $error;
 | 
|---|
 | 550 | }
 | 
|---|
 | 551 | 
 | 
|---|
 | 552 | sub ok {
 | 
|---|
 | 553 |     my $self = shift;
 | 
|---|
 | 554 |     $self->{"ok"}
 | 
|---|
 | 555 | }
 | 
|---|
 | 556 | 
 | 
|---|
 | 557 | sub error {
 | 
|---|
 | 558 |     my $self = shift;
 | 
|---|
 | 559 |     $self->{"error"};
 | 
|---|
 | 560 | }
 | 
|---|
 | 561 | 
 | 
|---|
 | 562 | sub exists {
 | 
|---|
 | 563 |     my $self = shift;
 | 
|---|
 | 564 |     $self->{"exists"}
 | 
|---|
 | 565 | }
 | 
|---|
 | 566 | 
 | 
|---|
 | 567 | sub input {
 | 
|---|
 | 568 |     my $self = shift;
 | 
|---|
 | 569 |     $self->{"qcinput"}
 | 
|---|
 | 570 | }
 | 
|---|
 | 571 | 
 | 
|---|
 | 572 | sub inputok {
 | 
|---|
 | 573 |     my $self = shift;
 | 
|---|
 | 574 |     $self->{"qcinput"}->ok();
 | 
|---|
 | 575 | }
 | 
|---|
 | 576 | 
 | 
|---|
 | 577 | sub energy {
 | 
|---|
 | 578 |     my $self = shift;
 | 
|---|
 | 579 |     $self->{"energy"}
 | 
|---|
 | 580 | }
 | 
|---|
 | 581 | 
 | 
|---|
 | 582 | sub s2norm {
 | 
|---|
 | 583 |     my $self = shift;
 | 
|---|
 | 584 |     $self->{"s2norm"}
 | 
|---|
 | 585 | }
 | 
|---|
 | 586 | 
 | 
|---|
 | 587 | sub s2matrix1norm {
 | 
|---|
 | 588 |     my $self = shift;
 | 
|---|
 | 589 |     $self->{"s2matrix1norm"}
 | 
|---|
 | 590 | }
 | 
|---|
 | 591 | 
 | 
|---|
 | 592 | sub degenerate {
 | 
|---|
 | 593 |     my $self = shift;
 | 
|---|
 | 594 |     $self->{"degenerate"}
 | 
|---|
 | 595 | }
 | 
|---|
 | 596 | 
 | 
|---|
 | 597 | sub d1mp2 {
 | 
|---|
 | 598 |     my $self = shift;
 | 
|---|
 | 599 |     $self->{"d1mp2"}
 | 
|---|
 | 600 | }
 | 
|---|
 | 601 | 
 | 
|---|
 | 602 | sub d2mp1 {
 | 
|---|
 | 603 |     my $self = shift;
 | 
|---|
 | 604 |     $self->{"d2mp1"}
 | 
|---|
 | 605 | }
 | 
|---|
 | 606 | 
 | 
|---|
 | 607 | sub s2matrixinfnorm {
 | 
|---|
 | 608 |     my $self = shift;
 | 
|---|
 | 609 |     $self->{"s2matrixinfnorm"}
 | 
|---|
 | 610 | }
 | 
|---|
 | 611 | 
 | 
|---|
 | 612 | sub t1norm {
 | 
|---|
 | 613 |     my $self = shift;
 | 
|---|
 | 614 |     $self->{"t1norm"}
 | 
|---|
 | 615 | }
 | 
|---|
 | 616 | 
 | 
|---|
 | 617 | sub t1matrix2norm {
 | 
|---|
 | 618 |     my $self = shift;
 | 
|---|
 | 619 |     $self->{"t1matrix2norm"}
 | 
|---|
 | 620 | }
 | 
|---|
 | 621 | 
 | 
|---|
 | 622 | sub optmolecule {
 | 
|---|
 | 623 |     my $self = shift;
 | 
|---|
 | 624 |     $self->{"optmolecule"};
 | 
|---|
 | 625 | }
 | 
|---|
 | 626 | 
 | 
|---|
 | 627 | sub gradient {
 | 
|---|
 | 628 |     my $self = shift;
 | 
|---|
 | 629 |     $self->{"grad"};
 | 
|---|
 | 630 | }
 | 
|---|
 | 631 | 
 | 
|---|
 | 632 | sub frequencies {
 | 
|---|
 | 633 |     my $self = shift;
 | 
|---|
 | 634 |     $self->{"freq"}
 | 
|---|
 | 635 | }
 | 
|---|
 | 636 | 
 | 
|---|
 | 637 | sub d1large_coef {
 | 
|---|
 | 638 |     my $self = shift;
 | 
|---|
 | 639 |     $self->{"d1large_coef"}
 | 
|---|
 | 640 | }
 | 
|---|
 | 641 | 
 | 
|---|
 | 642 | sub d1large_i {
 | 
|---|
 | 643 |     my $self = shift;
 | 
|---|
 | 644 |     $self->{"d1large_i"}
 | 
|---|
 | 645 | }
 | 
|---|
 | 646 | 
 | 
|---|
 | 647 | sub d1large_j {
 | 
|---|
 | 648 |     my $self = shift;
 | 
|---|
 | 649 |     $self->{"d1large_j"}
 | 
|---|
 | 650 | }
 | 
|---|
 | 651 | 
 | 
|---|
 | 652 | sub d1large_a {
 | 
|---|
 | 653 |     my $self = shift;
 | 
|---|
 | 654 |     $self->{"d1large_a"}
 | 
|---|
 | 655 | }
 | 
|---|
 | 656 | 
 | 
|---|
 | 657 | sub d1large_b {
 | 
|---|
 | 658 |     my $self = shift;
 | 
|---|
 | 659 |     $self->{"d1large_b"}
 | 
|---|
 | 660 | }
 | 
|---|
 | 661 | 
 | 
|---|
 | 662 | sub d1large_spin {
 | 
|---|
 | 663 |     my $self = shift;
 | 
|---|
 | 664 |     $self->{"d1large_spin"}
 | 
|---|
 | 665 | }
 | 
|---|
 | 666 | 
 | 
|---|
 | 667 | sub s2large_coef {
 | 
|---|
 | 668 |     my $self = shift;
 | 
|---|
 | 669 |     $self->{"s2large_coef"}
 | 
|---|
 | 670 | }
 | 
|---|
 | 671 | 
 | 
|---|
 | 672 | sub s2large_i {
 | 
|---|
 | 673 |     my $self = shift;
 | 
|---|
 | 674 |     $self->{"s2large_i"}
 | 
|---|
 | 675 | }
 | 
|---|
 | 676 | 
 | 
|---|
 | 677 | sub s2large_a {
 | 
|---|
 | 678 |     my $self = shift;
 | 
|---|
 | 679 |     $self->{"s2large_a"}
 | 
|---|
 | 680 | }
 | 
|---|
 | 681 | 
 | 
|---|
 | 682 | sub npacharge {
 | 
|---|
 | 683 |     my $self = shift;
 | 
|---|
 | 684 |     $self->{"npacharge"}
 | 
|---|
 | 685 | }
 | 
|---|
 | 686 | 
 | 
|---|
 | 687 | sub npashellpop {
 | 
|---|
 | 688 |     my $self = shift;
 | 
|---|
 | 689 |     $self->{"npashellpop"}
 | 
|---|
 | 690 | }
 | 
|---|
 | 691 | 
 | 
|---|
 | 692 | 1;
 | 
|---|