Hello, I have one important comment: The vectors coefficients in the nwchem output are incomplete! The default behaviour of nwchem is to print 10 first coefficients with value bigger than 0.15. For systems with many atoms it is not enough, usually its not even close. This behaviour is hard-coded in the nwchem source. To change this you must search each instance of movecs_print_anal in the source code and replace 0.15d0 for smaller value in appropriate calls. Furthermore you must change one loop in the src/ddscf/movecs_pr_anal.F file and around 200 line there will be loop: do klo = 0, min(n-1,9), 2 You must increase the range of this loop, for something more reasonable like: do klo = 0, min(n-1,199), 2 After recompiling the nwchem will print more coefficients and the gabedit will produce more reliable orbitals. Best regards, Karol Strutynski
So let's modify NWChem. I'll be modifying the 27th of June release of NWChem 6.1.1, which you'll obtain as Nwchem-6.1.1-src.2012-06-27.tar.gz from http://www.nwchem-sw.org/index.php/Download.
Change the number in red to something smaller (I tried 0.01d0) in the following files:
/src/ddscf/uhf.F
146 9611 continue
147 call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs,
148 $ 'UHF Final Alpha Molecular Orbital Analysis',
149 $ .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
150 $ .true., dbl_mb(k_occ))
151 call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs(2),
152 $ 'UHF Final Beta Molecular Orbital Analysis',
153 $ .true., dbl_mb(k_eval+nbf), oadapt, int_mb(k_irs+nmo),
154 $ .true., dbl_mb(k_occ+nbf)
/src/ddscf/scf_vec_guess.F
506 if (scftype.eq.'RHF' .or. scftype.eq.'ROHF') then
507 call movecs_print_anal(basis, 1,
508 & nprint, 0.15d0, g_movecs,
509 & 'ROHF Initial Molecular Orbital Analysis',
510 & .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
511 & .true., dbl_mb(k_occ))
512 else
513 nprint = min(nalpha+20,nmo)
514 call movecs_print_anal(basis, max(1,nbeta-20),
515 & nprint, 0.15d0, g_movecs,
516 & 'UHF Initial Alpha Molecular Orbital Analysis',
517 & .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
518 & .true., dbl_mb(k_occ))
519 call movecs_print_anal(basis, max(1,nbeta-20),
520 & nprint, 0.15d0, g_movecs(2),
521 & 'UHF Initial Beta Molecular Orbital Analysis',
522 & .true., dbl_mb(k_eval+nbf), oadapt, int_mb(k_irs+nmo),
523 & .true., dbl_mb(k_occ+nbf))
/src/ddscf/rohf.F
155 endif
156 call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs,
157 $ 'ROHF Final Molecular Orbital Analysis',
158 $ .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
159 $ .true., dbl_mb(k_occ))
/src/mcscf/mcscf.F
680 if (util_print('final vectors analysis', print_default))
681 $ call movecs_print_anal(basis,
682 $ max(1,nclosed-10), min(nbf,nclosed+nact+10),
683 $ 0.15d0, g_movecs, 'Analysis of MCSCF natural orbitals',
684 $ .true., dbl_mb(k_evals), .true., int_mb(k_sym),
685 $ .true., dbl_mb(k_occ))
686 c
/src/nwdft/scf_dft_cg/dft_cg_solve.F
166 call movecs_fix_phase(g_movecs(ispin))
167 call movecs_print_anal(basis, ilo, ihi, 0.15d0,
168 & g_movecs(ispin),blob,
169 & .true., dbl_mb(k_eval+(ispin-1)*nbf),
170 & oadapt, int_mb(k_irs+(ispin-1)*nbf),
171 & .true., dbl_mb(k_occ+(ispin-1)*nbf))
172 enddo
/src/nwdft/scf_dft/dft_scf.F
1736 call movecs_print_anal(ao_bas_han, ilo, ihi, 0.15d0,
1737 & g_movecs(ispin),
1738 & blob,
1739 & .true., dbl_mb(k_eval(ispin)), oadapt,
1740 & int_mb(k_ir+(ispin-1)*nbf_ao),
1741 & .true., dbl_mb(k_occ+(ispin-1)*nbf_ao))
/src/nwdft/scf_dft/dft_mxspin_ovlp.F
186 call movecs_print_anal(basis,int_mb(k_non),int_mb(k_non)
187 & ,0.15d0,g_alpha,'Alpha Orbitals without Beta Partners',
188 & .false., 0.0 ,.false., 0 , .false., 0 )
189 c
190 if (nct.GE.2) then
191 do i = 2,nct
192 ind = int_mb(k_non+i-1)
193 call movecs_print_anal(basis,ind,ind
194 & ,0.15d0,g_alpha,' ',
195 & .false., 0.0 ,.false., 0 , .false., 0 )
196 enddo
197 endif
352 c
353 call movecs_print_anal(basis, 1, nalp, 0.15d0, g_ualpha,
354 & 'Alpha Orb. w/o Beta Partners (after maxim. alpha/beta overlap)',
355 & .false., 0.0 ,.false., 0 , .false., 0 )
356 c
Otherwise once could presumably edit the header in ./src/ddscf/movecs_pr_anal.F directly and substitute thresh. At a minimum you should edit that file according to Karol's instructions: change the number in red below to e.g. 199./src/ddscf/movecs_pr_anal.F
198 do klo = 0, min(n-1,9), 2
199 khi = min(klo+1,n-1)
200 write(LuOut,2) (
201 $ int_mb(k_list+k)+1,
202 $ dbl_mb(k_vecs+int_mb(k_list+k)),
203 $ (byte_mb(k_tags+int_mb(k_list+k)*16+m),m=0,15),
204 $ k = klo,khi)
205 2 format(1x,2(i5,2x,f12.6,2x,16a1,4x))
206 enddo
Compilation
At this point you should be able to follow post 242. Briefly: Compiling NWChem 6.1.1 with Python on Debian Testing (Wheezy) and compile nwchem with python etc. Don't forget to edit /src/config/makefile.h for python support as shown in that post. Once you're done with that you can compare the GabEdit plots with and without the modification.
Alternatively, if you're simply making changes to a copy of nwchem that you've compiled before, you can speed thing up by a factor of ca 300 by following this post:
http://verahill.blogspot.com.au/2013/04/380-modifying-nwchem-code-without-full.html
The difference:
I ran a job on benzene as described in post 281. Visualising NWChem output with GabEdit. I chose to run use the ELF (electron localisation function) on output from the unmodified and modified nwchem binaries. It's a pretty big difference:
 |
| Original |
 |
| Modified |
