| 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;
 | 
|---|