Two bugs affecting ZMATRIX directive
Clicked A Few Times
5:57:01 AM PDT - Thu, May 12th 2016
Hi,
I think I've found two bugs in the routines dealing with internal coordinates. The bugs were identified in the latest svn r28223 of NWChem-6.6, but this seems the case also for NWChem-6.3 r24652.

  1. When using internal coordinates in the ZMATRIX directive, if any of the reference atoms is the first atom (atom index=1), and the 3rd coordinate is an angle (i.e., type BEND(2)) with orientation=+1, NWChem will return "SAME ATOM REFERRED TO MORE THAN ONCE ON THIS LINE". This is because the orientation of the last angle, i.e., the "+1", is perceived as an atom index.

  2. The values of BEND(2) coordinates are altered when the ZMatrix input is converted into Cartesian coordinates. The resulting geometry is distorted. This results in a distorted geometry. The problem is a minus sign missing in the code converting internal to Cartesian coords.


None of these bugs can be detected running the test QA/tests/geom_zmatrix/geom_zmatrix.nw. In fact, the only BEND(2) angle in this test has orientation -1, thus escapes bug 1, and its value does not introduce any significant out-of-plane angle, thus the geometrical distortion from bug 2 is negligible.

The following geometry could be more effective for testing the routines against these two bugs.
GEOMETRY
  ZMATRIX
    C
    O  1 dst1
    C  1 dst2 2 ang1
    H  1 dst3 2 ang2 3 ang3 1
    H  1 dst4 2 ang4 3 ang5 -1
    H  2 dst5 1 ang6 3 tor1
    H  3 dst6 1 ang7 2 tor2
    H  3 dst7 1 ang8 7 ang9 1
    H  3 dst8 1 ang10 7 ang11 -1
    VARIABLES
     ang10  109.47042118
     dst8  1.08998960
     ang5  109.47042118
     ang4  109.46955544
     ang11  109.46955544
     dst6  1.09000000
     ang3  109.47042118
     ang8  109.47042118
     ang2  109.46955544
     dst5  0.98005555
     dst3  1.08998960
     dst1  1.42000000
     ang6  109.47213930
     ang9  109.46955544
     dst2  1.53000204
     ang7  109.47119359
     dst4  1.08998960
     dst7  1.08998960
     ang1  109.47119359
     tor2  180.00000000
     tor1  180.00000000
  END
END


Finally, here is the patch to fix the two bugs. It also removes the case sensitivity towards the END/ZEND label of the ZMATRIX directive. This is to make it consistent with the general syntax for directives. 

Index: geom/geom_hnd.F
===================================================================
--- geom/geom_hnd.F	(revision 28223)
+++ geom/geom_hnd.F	(working copy)
@@ -3918,7 +3918,7 @@
          WRITE(IW,9995) IAT,(ZMT(I,IAT),I=1,5),(  ZLST(I,IAT),I=1,3),
      1                                         (FRZVAL(I,IAT),I=1,3)
       ENDIF
-      DO K=1,5
+      DO K=1,4
         IF (ZMT(K,IAT).NE.0) THEN
           DO J=1,K-1
             IF (ZMT(K,IAT).EQ.ZMT(J,IAT)) THEN
@@ -4250,7 +4250,7 @@
                      numd=(COS(BET)- COS(ALP)* COS(GAM))/
      /                    (SIN(ALP)* SIN(GAM))
                      if(numd.gt.1d0) numd=1d0
-                     if(numd.lt.1d0) numd=-1d0
+                     if(numd.lt.-1d0) numd=-1d0
                        THETA=ACOS(numd)
                         IF(NZMT(5,IAT).EQ.-1) THEN
                            THETA=-THETA
Index: geom/geom_input.F
===================================================================
--- geom/geom_input.F	(revision 28223)
+++ geom/geom_input.F	(working copy)
@@ -3446,8 +3446,10 @@
          IF(NUMWRD.GT.1) THEN
             GO TO 110
          ELSE
-            IF((NUMCHR(1).EQ.6.AND.PRSWRD(1)(2:5).EQ.CHREND(1:4)).OR.
-     1         (NUMCHR(1).EQ.5.AND.PRSWRD(1)(2:4).EQ.CHREND(2:4))) THEN
+            IF((NUMCHR(1).EQ.6.AND.
+     &         INP_COMPARE(.FALSE.,CHREND(1:4),PRSWRD(1)(2:5))).OR.
+     &         (NUMCHR(1).EQ.5.AND.
+     &         INP_COMPARE(.FALSE.,CHREND(2:4),PRSWRD(1)(2:4)))) THEN
                NAT = IAT -1
                GO TO 140
             ELSE
@@ -3491,8 +3493,10 @@
             ENDIF
             GO TO 120
          ELSE
-            IF((NUMCHR(1).EQ.6.AND.PRSWRD(1)(2:5).EQ.CHREND(1:4)).OR.
-     1         (NUMCHR(1).EQ.5.AND.PRSWRD(1)(2:4).EQ.CHREND(2:4))) THEN
+            IF((NUMCHR(1).EQ.6.AND.
+     &         INP_COMPARE(.FALSE.,CHREND(1:4),PRSWRD(1)(2:5))).OR.
+     &         (NUMCHR(1).EQ.5.AND.
+     &         INP_COMPARE(.FALSE.,CHREND(2:4),PRSWRD(1)(2:4)))) THEN
                NVAR = IVAR -1
                GO TO 140
             ELSE
@@ -3526,8 +3530,10 @@
                FRZVAR(  IVAR)=.TRUE.
             GO TO 130
          ELSEIF(NUMWRD.EQ.1) THEN
-            IF((NUMCHR(1).EQ.6.AND.PRSWRD(1)(2:5).EQ.CHREND(1:4)).OR.
-     1         (NUMCHR(1).EQ.5.AND.PRSWRD(1)(2:4).EQ.CHREND(2:4))) THEN
+            IF((NUMCHR(1).EQ.6.AND.
+     &         INP_COMPARE(.FALSE.,CHREND(1:4),PRSWRD(1)(2:5))).OR.
+     &         (NUMCHR(1).EQ.5.AND.
+     &         INP_COMPARE(.FALSE.,CHREND(2:4),PRSWRD(1)(2:4)))) THEN
             ELSE
                WRITE(IW,8880)
                CALL HND_HNDERR(3,ERRMSG)


Now we can have fun with internal coordinates!

Forum Vet
7:02:54 PM PDT - Thu, May 12th 2016
Thanks
Thanks for providing this patch.
I have checked in to the svn repository.
Could you provide your name and affiliation for proper attribution of the contribution?

Cheers, Edo

Clicked A Few Times
1:36:10 AM PDT - Fri, May 13th 2016
Sure! I'm Marco Foscato, University of Bergen


Forum >> NWChem's corner >> General Topics