12 February 2013

337. Modifying Nwchem 6.1.1 to work with GabEdit

Karol Strutynski left the following comment on a post about NWChem and Gabedit:

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

4 comments:

  1. Hello:

    Thank you so much for this valuable information.

    But under the heading of "Compilation" above, the directions lead to a pointer to compile the entire NWChem source. But only one or two subroutines are modified. Shouldn't there be a more surgical way of proceeding after making the suggested changes? Actually this would make a great separate post. That is, provide a procedure to recompile and link when small changes to NWChem are made. I've done this, but I can't remember.

    Kind regards,

    Angelo

    ReplyDelete
    Replies
    1. Angelo,
      I seem to remember seeing Bert describe how to do it, but also don't remember. I'll work it out and make it a separate post (might take a few days). Thanks for the suggestion!

      Delete
    2. Havne't written it up yet, but it turned out to be easy (looked at this page: http://xray.isc.kharkov.com/ext_docs/NWChem/prog/node12.html)

      1. set all the environment variables (e.g. NWCHEM_TOP, FC, etc.) as done during a full build
      2. Go to
      a) src/ddscf, and run make
      b) src/nwdft/scf_dft and run make
      c) src/mcscf and run make
      3. Go to src/ and run make link

      Total time is 0.3+2+0.8+3=6.1 seconds, instead of ca 30 minutes for a full recompile. I've tested it and have confirmed that it works.

      Delete
    3. And it's here now: http://verahill.blogspot.com.au/2013/04/380-modifying-nwchem-code-without-full.html

      Delete