Wilson-Amos-Handy method for calculation of NMR shifts

Gets Around

5:15:21 PM PST - Thu, Nov 8th 2012
Dear all, My, probably naive, question is this: is it possible to use NWChem in its current state for the Wilson-Amos-Handy method to calculating NMR shifts? The approach is described in Chem. Phys. Letters 1999, 312, 475-484 and J. Chem. Phys, 2000, 113(8), 2983-2989 and 'excludes the exact exchange from the calculation of the first-order perturbed wave function so that the coupled-perturbed KS can be diagonalized and doesn't need to be iteratively solved' (paraphrasing Magyarfalvi and Pulay's description on page 1351 in J. Chem. Phys. 2003, 119, 1350-1356). At a first glance it may look like the method only involves a slight re-parametrisation (xc HFexch 0.05 slater 0.95 becke88 nonlocal 0.72 perdew91 0.81 wvn_5 1), but the functional alone performs as poorly as b3lyp, acm etc. and there are more significant changes in the approach. To my knowledge this has only been implemented in PQS so far. I'm interested in the method given the reported excellent accuracy in the prediction of the 17O NMR shifts ('simple' test cases like CO and H2CO are not handled well by giao/b3lyp) for at least first and second row elements. Is it possible to reproduce the Wilson-Amos-Handy method in the current version of NWChem, or would it involve non-trivial changes/additions to the code? Best regards, Andy EDIT: fixed the sloppily type exchange correlation functional definition

Forum Vet

6:08:51 PM PST - Thu, Nov 8th 2012
If I read the papers correctly, for B3LYP the functional woudl have the following form: XC HFexch 0.05 slater 0.95 vwn_1_rpa 0.19 lyp 0.81 becke88 nonlocal 0.72 Explicitly not including the HFexch in the NMR response calculation would require some minor modifications to make it feasible. Bert Quote:Ohlincha Nov 9th 12:15 am Dear all, My, probably naive, question is this: is it possible to use NWChem in its current state for the Wilson-Amos-Handy method to calculating NMR shifts? The approach is described in Chem. Phys. Letters 1999, 312, 475-484 and J. Chem. Phys, 2000, 113(8), 2983-2989 and 'excludes the exact exchange from the calculation of the first-order perturbed wave function so that the coupled-perturbed KS can be diagonalized and doesn't need to be iteratively solved' (paraphrasing Magyarfalvi and Pulay's description on page 1351 in J. Chem. Phys. 2003, 119, 1350-1356). At a first glance it may look like the method only involves a slight re-parametrisation (xc HFexch 0.05 slater 0.95 becke88 0.72 perdew 91 wvn_5 1), but the functional alone performs as poorly as b3lyp, acm etc. and there are more significant changes in the approach. To my knowledge this has only been implemented in PQS so far. I'm interested in the method given the reported excellent accuracy in the prediction of the 17O NMR shifts ('simple' test cases like CO and H2CO are not handled well by giao/b3lyp) for at least first and second row elements. Is it possible to reproduce the Wilson-Amos-Handy method in the current version of NWChem, or would it involve non-trivial changes/additions to the code? Best regards, Andy

Gets Around

11:03:18 PM PST - Thu, Nov 8th 2012
Bert, The PQS energies and the remarks of Pulay (which is probably a bit more relevant) point towards the original authors using the form defined by Becke in J . Chem. Phys. 1993, 98(7), 5648-5652 (i.e. what's included in NWChem as the 'acm' functional) instead of the Gaussian '92 version (what's included as b3lyp in NWChem and most other programmes), so they most likely used VWN_5 and PW91 instead of VWN_1_RPA and (LYP-VWN_1_RPA). Anyway, thank you for your prompt response. Perhaps this is something that may be incorporated somewhere down the line? /Andy

Forum Vet

10:55:44 AM PST - Fri, Nov 9th 2012
We'll put adding a keyword and functionality to exclude exact exchange from the response calculation on our list. Thanks, Bert Quote:Ohlincha Nov 9th 6:03 am Bert, The PQS energies and the remarks of Pulay (which is probably a bit more relevant) point towards the original authors using the form defined by Becke in J . Chem. Phys. 1993, 98(7), 5648-5652 (i.e. what's included in NWChem as the 'acm' functional) instead of the Gaussian '92 version (what's included as b3lyp in NWChem and most other programmes), so they most likely used VWN_5 and PW91 instead of VWN_1_RPA and (LYP-VWN_1_RPA). Anyway, thank you for your prompt response. Perhaps this is something that may be incorporated somewhere down the line? /Andy

Gets Around

6:38:27 PM PST - Wed, Nov 14th 2012
Cheers Bert. Let me know if/when it makes it into the development version and I'll be happy to test it against the published data set. Hopefully such a switch will be enough to make it work - I tried editing (in NWChem 6.1) a number of files in src/property: aoresponse_giao_rhs.F, giao_b1_movecs.F, giaofock.F, hnd_giaox.F, hnd_giaox_zora.F, int_giao_1ega.F and changing to xfac = 0.0d0 c if (use_theory.eq.'dft') xfac = bgj_kfac() but then again I would claim to know what I'm doing and it obviously still iteratively solves the matrix rather than diagonalise it. /Andy

Forum Vet

9:45:38 AM PST - Thu, Nov 15th 2012
Yes, it would still do it iteratively, but using only one iteration as the linear equation would be now uncoupled. Effectively, it diagonalizes it in one iteration. Hence, if that's what you see when you try it, you may have created a solution for yourself. Bert Quote:Ohlincha Nov 15th 1:38 am Cheers Bert. Let me know if/when it makes it into the development version and I'll be happy to test it against the published data set. Hopefully such a switch will be enough to make it work - I tried editing (in NWChem 6.1) a number of files in src/property: aoresponse_giao_rhs.F, giao_b1_movecs.F, giaofock.F, hnd_giaox.F, hnd_giaox_zora.F, int_giao_1ega.F and changing to xfac = 0.0d0 c if (use_theory.eq.'dft') xfac = bgj_kfac() but then again I would claim to know what I'm doing and it obviously still iteratively solves the matrix rather than diagonalise it. /Andy

Gets Around

10:07:28 PM PST - Thu, Nov 15th 2012

Bert,
I don't see much difference between the original 6.1 and my modified 6.1 code.

The number of iterations are the same, and while there's a small change in the computed shift, it's small relative what should be obtained if the WIlson-Amos-Handy method has been successfully implemented (isotropic shifts -- C: 4.7 ppm and O: -44 ppm; see first post for references)

Here's my test case

title "co nmr"
echo

geometry
c   0   0   0
o   0   0   1.13
end

basis 
    * library "dzp (dunning)"
end

property
    shielding
end

dft
 direct
 grid fine
 mult 1
 xc HFexch 0.05 slater 0.95 becke88 nonlocal 0.72 vwn_5 1 perdew91 0.81
end

task dft property

Running it using my modified binary (v6.1) I get:

600    iter   nsub   residual    time
601    ----  ------  --------  ---------
602      1      3    1.70D-01       4.9
603      2      6    1.70D-03       5.6
604      3      9    2.93D-05       6.3

For C I get

621            isotropic =      -2.8388$
622           anisotropy =     413.3050$

and for O I get

650            isotropic =     -64.1336$
651           anisotropy =     711.1241$

and with the original binary I see

600    iter   nsub   residual    time$
601    ----  ------  --------  ---------$
602      1      3    1.69D-01       5.0$
603      2      6    1.84D-03       5.8$
604      3      9    2.47D-05       6.5$

For C I get

621            isotropic =       5.1483$
622           anisotropy =     401.3242$

and for O I get

650            isotropic =     -60.6964$
651           anisotropy =     705.9682$

when I do

cat *|grep "xfac ="

I can't find any unedited examples -- i.e. I've already changed them all to 0 and have commented out "if (use_theory.eq.'dft') xfac = bgj_kfac()"

I guess the most plausible explanation would be that there's more to these NMR calculations than it at first appears to be -- but then again, as you point out there should be a change in the number of iterations and there isn't.

Forum Vet

3:40:13 PM PST - Sat, Nov 17th 2012

Yeah, because you're only taking the HF exchange out of the RHS, it's still included in the CPHF calculations later on. Try the attached hnd_giaox.F. I'm doing this while on travel. I think this may work, but I have not tested it.

Bert

[QUOTE=Ohlincha Nov 16th 5:07 am]Bert,

I don't see much difference between the original 6.1 and my modified 6.1 code.

The number of iterations are the same, and while there's a small change in the computed shift, it's small relative what should be obtained if the WIlson-Amos-Handy method has been successfully implemented (isotropic shifts -- C: 4.7 ppm and O: -44 ppm; see first post for references)

Here's my test case

title "co nmr"
echo

geometry
c   0   0   0
o   0   0   1.13
end

basis 
    * library "dzp (dunning)"
end

property
    shielding
end

dft
 direct
 grid fine
 mult 1
 xc HFexch 0.05 slater 0.95 becke88 nonlocal 0.72 vwn_5 1 perdew91 0.81
end

task dft property

Running it using my modified binary (v6.1) I get:

600    iter   nsub   residual    time
601    ----  ------  --------  ---------
602      1      3    1.70D-01       4.9
603      2      6    1.70D-03       5.6
604      3      9    2.93D-05       6.3

For C I get

621            isotropic =      -2.8388$
622           anisotropy =     413.3050$

and for O I get

650            isotropic =     -64.1336$
651           anisotropy =     711.1241$

and with the original binary I see

600    iter   nsub   residual    time$
601    ----  ------  --------  ---------$
602      1      3    1.69D-01       5.0$
603      2      6    1.84D-03       5.8$
604      3      9    2.47D-05       6.5$

For C I get

621            isotropic =       5.1483$
622           anisotropy =     401.3242$

and for O I get

650            isotropic =     -60.6964$
651           anisotropy =     705.9682$

when I do

cat *|grep "xfac ="

Gets Around

2:58:49 PM PST - Mon, Nov 19th 2012

Bert,
I couldn't kompile with the file you provided (debian, gcc 4.7),but with the following modifications it worked fine:

16c16
< c#include "../nwdft/include/cdft.fh"
---
> #include "../nwdft/include/cdft.fh"
33c33
<       double precision jfac(3),kfac(3),a(6),axs(3,3),eig(3),xfac(1)
---
>       double precision jfac(3),kfac(3),a(6),axs(3,3),eig(3)
42d41
<       logical do_wah

To my nwchem input I added

set dft:wah true

which I guess should be right:

cat co.db|strings|grep dft
   dft:wah
   !rtdb!dft:wah

but the shifts are the same as for the unmodified code (see below), and the number of iterations are not affected.

This is the case both when I modify xfac in the other files in src/property, and when I don't.

I also tried forcing wah by removing the if construct (i.e. same as if do_wah is true every time)

16c16
< c#include "../nwdft/include/cdft.fh"
---
> #include "../nwdft/include/cdft.fh"
33c33
<       double precision jfac(3),kfac(3),a(6),axs(3,3),eig(3),xfac(1)
---
>       double precision jfac(3),kfac(3),a(6),axs(3,3),eig(3)
76,78c76,78
< c      if (.not. rtdb_get(rtdb, 'dft:wah', mt_log, 1, do_wah))
< c     $   do_wah = .false.
< c      if (do_wah) then
---
>       if (.not. rtdb_get(rtdb, 'dft:wah', mt_log, 1, do_wah))
>      $   do_wah = .false.
>       if (do_wah) then
81c81
< c      endif
---
>       endif
872c872
< c      if (do_wah) then
---
>       if (do_wah) then
874c874
< c      endif
---
>       endif

But it didn't make any difference.

For either version of the code I get

For Carbon:

           isotropic =       5.1483
          anisotropy =     401.3242

For Oxygen:

           isotropic =     -60.6964
          anisotropy =     705.9682

I've tried a number of permutations but with little luck in terms of isotopic shift changes.

PS It's also become more difficult to post on the forum in the past week or two owing to the 'The specified URL cannot be found" issue.

Forum Vet

10:13:05 AM PST - Wed, Nov 21st 2012

Your fixing is never ever going to work. Key is the modification of xfac in the data block that is in cdft.fh . This is what is used by the CPHF routine and outside of hnd_giaox.F. You cannot simply fix this within this one routine. CPHF does the response calculation using the RHS constructed in hnd_giaox.F. There it also uses the HF Exchange.

I can see why it may not have compiled, forgot to define do_wah as a logical.

Hence, in the hnd_giaox.F file I provided, add somewhere at the top where all the variables are defined:

     logical  do_wah

Also, in the input add the following (as you already figured out):

     set dft:wah  .true.

If it doesn't compile, get me the error message.

Bert

Quote:Bert Nov 17th 10:40 pm

[QUOTE=Ohlincha Nov 16th 5:07 am]Bert,

title "co nmr"
echo

geometry
c   0   0   0
o   0   0   1.13
end

basis 
    * library "dzp (dunning)"
end

property
    shielding
end

dft
 direct
 grid fine
 mult 1
 xc HFexch 0.05 slater 0.95 becke88 nonlocal 0.72 vwn_5 1 perdew91 0.81
end

task dft property

Running it using my modified binary (v6.1) I get:

600    iter   nsub   residual    time
601    ----  ------  --------  ---------
602      1      3    1.70D-01       4.9
603      2      6    1.70D-03       5.6
604      3      9    2.93D-05       6.3

For C I get

621            isotropic =      -2.8388$
622           anisotropy =     413.3050$

and for O I get

650            isotropic =     -64.1336$
651           anisotropy =     711.1241$

and with the original binary I see

600    iter   nsub   residual    time$
601    ----  ------  --------  ---------$
602      1      3    1.69D-01       5.0$
603      2      6    1.84D-03       5.8$
604      3      9    2.47D-05       6.5$

For C I get

621            isotropic =       5.1483$
622           anisotropy =     401.3242$

and for O I get

650            isotropic =     -60.6964$
651           anisotropy =     705.9682$

when I do

cat *|grep "xfac ="

Gets Around

2:58:48 PM PST - Wed, Nov 21st 2012

giaofock.F:371:0: warning: ‘ntot’ may be used uninitialized in this function [-Wmaybe-uninitialized]
giaofock.F:508:0: warning: ‘factor_cam’ may be used uninitialized in this function [-Wmaybe-uninitialized]
../nwdft/include/cdft.fh:82.29:
    Included at hnd_giaox.F:17:

     &     nrinc, nrmax, geom, ncenters, nbf, nradpts, nang_leb_pts,
                               1
Error: COMMON attribute conflicts with DUMMY attribute in 'geom' at (1)
hnd_giaox.F:20.18:

      integer geom    ! [input] geometry handle                         
                  1
Error: Symbol 'geom' at (1) already has basic type of INTEGER
hnd_giaox.F:22.56:

      integer nclosed(2), nopen(2), nvirt(2), ndens, nbf, nmo           
                                                        1
Error: Symbol 'nbf' at (1) already has basic type of INTEGER
hnd_giaox.F:45.23:

      logical     oskel, status, debug                                  
                       1
Error: Symbol 'oskel' at (1) already has basic type of LOGICAL
hnd_giaox.F:81.13:

         xfac(1) = 0.0d0                                                
             1
Error: 'xfac' at (1) is not a variable
hnd_giaox.F:874.13:

         xfac(1) = kfac_org                                             
             1
Error: 'xfac' at (1) is not a variable
hnd_giaox.F:66.11:

      debug = .false. .and. (ga_nodeid().eq.0) ! special debugging      
           1
Error: Symbol debug at (1) has no IMPLICIT type
hnd_giaox.F:144.56:

      call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp,        
                                                        1
Error: Symbol ndens at (1) has no IMPLICIT type
hnd_giaox.F:140.46:

     &                      dbl_mb(k_eval),nmo)                         
                                              1
Error: Symbol nmo at (1) has no IMPLICIT type
hnd_giaox.F:138.61:

      call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen,          
                                                             1
Error: Symbol nopen at (1) has no IMPLICIT type
hnd_giaox.F:101.12:

      status = rtdb_parallel(.true.)                                    
            1
Error: Symbol status at (1) has no IMPLICIT type
hnd_giaox.F:138.55:

      call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen,          
                                                       1
Error: Symbol nclosed at (1) has no IMPLICIT type
hnd_giaox.F:161.15:

      ahi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:170.47:

      if(.not.ga_create(MT_DBL,nclosed(1)*nvirt(1),3,RHS,-1,-1,g_rhs))
                                               1
Error: Function nvirt at (1) has no IMPLICIT type
hnd_giaox.F:180.31:

      bhi(1) = nclosed(1)*nvirt(1)                                      
                               1
Error: Function nvirt at (1) has no IMPLICIT type
hnd_giaox.F:202.15:

      alo(1) = nclosed(1)+1                                             
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:206.18:

      do iocc = 1, nclosed(1)                                           
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:220.15:

      ahi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:230.15:

      ahi(1) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:264.15:

      dhi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:278.18:

         ahi(2) = nclosed(1)                                            
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:283.18:

         ahi(1) = nclosed(1)                                            
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:284.18:

         bhi(2) = nclosed(1)                                            
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:328.15:

      ahi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:340.15:

      clo(2) = nclosed(1)+1                                             
               1
Error: Function nclosed at (1) has no IMPLICIT type
Fatal Error: Error count reached limit of 25.

Im building the development version from October using the following script on debian (gcc 4.7):

export LARGE_FILES=TRUE
export TCGRSH=/usr/bin/ssh
export NWCHEM_TOP=`pwd`
export NWCHEM_TARGET=LINUX64
export NWCHEM_MODULES="all"
export BLASOPT="-L/opt/openblas/lib -lopenblas"
export USE_MPI=y
export USE_MPIF=y
export USE_MPIF4=y
export MPI_LOC=/usr/lib/openmpi/lib
export MPI_INCLUDE=/usr/lib/openmpi/include
export LIBRARY_PATH=$LIBRARY_PATH:/usr/lib/openmpi/lib:/opt/openblas/lib
export LIBMPI="-lmpi -lopen-rte -lopen-pal -ldl -lmpi_f77 -lpthread"
cd $NWCHEM_TOP/src
export FC=gfortran
make clean
make nwchem_config
make FC=gfortran 1> make.log 2>make.err
cd ../contrib
./getmem.nwchem

It builds the unmodified code without any issues.

(The "URL cannot be found" is apparently also triggered by too many apostrofes in a code bracket environment)

Forum Vet

5:55:52 PM PST - Wed, Nov 21st 2012

OK. As I said, I gave you a possible solution but I had no chance to compile it remotely. Back home now, I made some changes. Try the attached hnd_giaox.F.

Bert

Quote:Ohlincha Nov 21st 9:58 pm

giaofock.F:371:0: warning: ‘ntot’ may be used uninitialized in this function [-Wmaybe-uninitialized]
giaofock.F:508:0: warning: ‘factor_cam’ may be used uninitialized in this function [-Wmaybe-uninitialized]
../nwdft/include/cdft.fh:82.29:
    Included at hnd_giaox.F:17:

     &     nrinc, nrmax, geom, ncenters, nbf, nradpts, nang_leb_pts,
                               1
Error: COMMON attribute conflicts with DUMMY attribute in 'geom' at (1)
hnd_giaox.F:20.18:

      integer geom    ! [input] geometry handle                         
                  1
Error: Symbol 'geom' at (1) already has basic type of INTEGER
hnd_giaox.F:22.56:

      integer nclosed(2), nopen(2), nvirt(2), ndens, nbf, nmo           
                                                        1
Error: Symbol 'nbf' at (1) already has basic type of INTEGER
hnd_giaox.F:45.23:

      logical     oskel, status, debug                                  
                       1
Error: Symbol 'oskel' at (1) already has basic type of LOGICAL
hnd_giaox.F:81.13:

         xfac(1) = 0.0d0                                                
             1
Error: 'xfac' at (1) is not a variable
hnd_giaox.F:874.13:

         xfac(1) = kfac_org                                             
             1
Error: 'xfac' at (1) is not a variable
hnd_giaox.F:66.11:

      debug = .false. .and. (ga_nodeid().eq.0) ! special debugging      
           1
Error: Symbol debug at (1) has no IMPLICIT type
hnd_giaox.F:144.56:

      call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp,        
                                                        1
Error: Symbol ndens at (1) has no IMPLICIT type
hnd_giaox.F:140.46:

     &                      dbl_mb(k_eval),nmo)                         
                                              1
Error: Symbol nmo at (1) has no IMPLICIT type
hnd_giaox.F:138.61:

      call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen,          
                                                             1
Error: Symbol nopen at (1) has no IMPLICIT type
hnd_giaox.F:101.12:

      status = rtdb_parallel(.true.)                                    
            1
Error: Symbol status at (1) has no IMPLICIT type
hnd_giaox.F:138.55:

      call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen,          
                                                       1
Error: Symbol nclosed at (1) has no IMPLICIT type
hnd_giaox.F:161.15:

      ahi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:170.47:

      if(.not.ga_create(MT_DBL,nclosed(1)*nvirt(1),3,RHS,-1,-1,g_rhs))
                                               1
Error: Function nvirt at (1) has no IMPLICIT type
hnd_giaox.F:180.31:

      bhi(1) = nclosed(1)*nvirt(1)                                      
                               1
Error: Function nvirt at (1) has no IMPLICIT type
hnd_giaox.F:202.15:

      alo(1) = nclosed(1)+1                                             
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:206.18:

      do iocc = 1, nclosed(1)                                           
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:220.15:

      ahi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:230.15:

      ahi(1) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:264.15:

      dhi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:278.18:

         ahi(2) = nclosed(1)                                            
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:283.18:

         ahi(1) = nclosed(1)                                            
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:284.18:

         bhi(2) = nclosed(1)                                            
                  1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:328.15:

      ahi(2) = nclosed(1)                                               
               1
Error: Function nclosed at (1) has no IMPLICIT type
hnd_giaox.F:340.15:

      clo(2) = nclosed(1)+1                                             
               1
Error: Function nclosed at (1) has no IMPLICIT type
Fatal Error: Error count reached limit of 25.

Im building the development version from October using the following script on debian (gcc 4.7):

export LARGE_FILES=TRUE
export TCGRSH=/usr/bin/ssh
export NWCHEM_TOP=`pwd`
export NWCHEM_TARGET=LINUX64
export NWCHEM_MODULES="all"
export BLASOPT="-L/opt/openblas/lib -lopenblas"
export USE_MPI=y
export USE_MPIF=y
export USE_MPIF4=y
export MPI_LOC=/usr/lib/openmpi/lib
export MPI_INCLUDE=/usr/lib/openmpi/include
export LIBRARY_PATH=$LIBRARY_PATH:/usr/lib/openmpi/lib:/opt/openblas/lib
export LIBMPI="-lmpi -lopen-rte -lopen-pal -ldl -lmpi_f77 -lpthread"
cd $NWCHEM_TOP/src
export FC=gfortran
make clean
make nwchem_config
make FC=gfortran 1> make.log 2>make.err
cd ../contrib
./getmem.nwchem

It builds the unmodified code without any issues.

(The "URL cannot be found" is apparently also triggered by too many apostrofes in a code bracket environment)

Gets Around

9:58:04 PM PST - Wed, Nov 21st 2012

Bert,
Thank you for putting all this effort into this.

NWChem now compiles fine and without any other modification.

Running a normal calculation works fine and gives the expected results.

Running a calculation with

set dft:wah .true.

however gives a fatak error

                               NWChem CPHF Module
                                ------------------


  scftype          =     RHF 
  nclosed          =        7
  nopen            =        0
  variables        =      175
  # of vectors     =        3
  tolerance        = 0.10D-03
  level shift      = 0.00D+00
  max iterations   =       50
  max subspace     =       30

 SCF residual:   0.36977616114820894     
 ------------------------------------------------------------------------
 cphf_solve2:SCF residual greater than 1d-2        0
 ------------------------------------------------------------------------
 This error has not yet been assigned to a category
 ------------------------------------------------------------------------
 For more information see the NWChem manual at http://nwchemgit.github.io/index.php/NWChem_Documentation


 For further details see manual section: No section for this category                                                                                                                                                                                                                                   
0:0:cphf_solve2:SCF residual greater than 1d-2:: 0
(rank:0 hostname:beryllium pid:1125):ARMCI DASSERT fail. ../../ga-5-1/armci/src/common/armci.c:ARMCI_Error():208 cond:0

I presume that the residual will somehow corresponds to the HF exchange and that there's an additional modification which is needed somewhere.

Forum Vet

10:29:30 AM PST - Mon, Nov 26th 2012

Searching for the error message you would have been able to find where this is happening. Try fixing it in cphf/cphf_solve2.F. Don't know if it will work as it should but worth a short. After that, I have to give up on finding quick solution.

Bert

Quote:Ohlincha Nov 22nd 4:58 am

set dft:wah .true.

however gives a fatak error

                               NWChem CPHF Module
                                ------------------


  scftype          =     RHF 
  nclosed          =        7
  nopen            =        0
  variables        =      175
  # of vectors     =        3
  tolerance        = 0.10D-03
  level shift      = 0.00D+00
  max iterations   =       50
  max subspace     =       30

 SCF residual:   0.36977616114820894     
 ------------------------------------------------------------------------
 cphf_solve2:SCF residual greater than 1d-2        0
 ------------------------------------------------------------------------
 This error has not yet been assigned to a category
 ------------------------------------------------------------------------
 For more information see the NWChem manual at http://nwchemgit.github.io/index.php/NWChem_Documentation


 For further details see manual section: No section for this category                                                                                                                                                                                                                                   
0:0:cphf_solve2:SCF residual greater than 1d-2:: 0
(rank:0 hostname:beryllium pid:1125):ARMCI DASSERT fail. ../../ga-5-1/armci/src/common/armci.c:ARMCI_Error():208 cond:0

I presume that the residual will somehow corresponds to the HF exchange and that there's an additional modification which is needed somewhere.

Gets Around

8:48:12 PM PST - Mon, Nov 26th 2012

Bert,
I can get it to run by simply suppressing the error as shown below, but that didn't solve the issue. I guess it's time to sit down and figure out what's executed and in what order.

52d51
<       logical do_wah
58,60d56
< 
< 
< 
244,248d239
<        
< c     check for Wilson-Amos-Handy NMR computation request
<       if (.not. rtdb_get(rtdb, 'dft:wah', mt_log, 1, do_wah)) 
<      $  do_wah=.false.
<       if (.not. (do_wah)) then
250,251c241,242
<           call ga_sync()
<           call errquit('cphf_solve2:SCF residual greater than 1d-2',
---
>          call ga_sync()
>          call errquit('cphf_solve2:SCF residual greater than 1d-2',
253,254c244
<        endif
<        endif
---
>       endif

Execution of the carbon monoxide job results in

 SCF residual:   0.36977288699976185     


Iterative solution of linear equations
  No. of variables      175
  No. of equations        3
  Maximum subspace       30
        Iterations       50
       Convergence  1.0D-04
        Start time      3.4


   iter   nsub   residual    time
   ----  ------  --------  ---------
     1      3    4.20D-02       4.1
     2      6    7.11D-04       4.8
     3      9    1.72D-05       5.5

which gives for O


           isotropic =      -2.5674
          anisotropy =     412.8978

and for C

           isotropic =     -63.0277
          anisotropy =     709.465

Forum >> NWChem's corner >> NWChem functionality