Write out total forces to out.xml + Bugfixes

1. The total forces are now written out to out.xml iff l_f=T.
2. There were bugs in some array dimensionings in force_a21. Corrected.
3. Some cleanup in cdnval, force_a21, and force_a21_lo.

cdn/cdnval.F90

@@ -782,7 +782,7 @@ CONTAINS
                                we,ispin,noccbd,usdus,acof(:,0:,:,ispin),&
                                bcof(:,0:,:,ispin),e1cof,e2cof, acoflo,bcoflo, results,f_a12)
 #endif
-                CALL force_a21(atoms,dimension,noccbd,sym,atoms%nlod*(atoms%nlod+1)/2,&
+                CALL force_a21(atoms,dimension,noccbd,sym,&
                                oneD,cell,we,ispin,epar(0:,:,ispin),noccbd,eig,usdus,acof(:,0:,:,ispin),&
                                bcof(:,0:,:,ispin),ccof(-atoms%llod:,:,:,:,ispin), aveccof,bveccof,cveccof,&
                                results,f_a21,f_b4)

force/force_a21.F90

 MODULE m_forcea21
 CONTAINS
   SUBROUTINE force_a21(&
-       atoms,dimension,nobd,sym, loplod,oneD,cell,&
+       atoms,dimension,nobd,sym,oneD,cell,&
        we,jsp,epar,ne,eig,usdus,&
        acof,bcof,ccof,aveccof,bveccof,cveccof, results,f_a21,f_b4)
 
@@ -38,8 +38,7 @@ CONTAINS
     TYPE(t_usdus),INTENT(IN)     :: usdus
     !     ..
     !     .. Scalar Arguments ..
-    INTEGER, INTENT (IN) :: nobd  
-    INTEGER, INTENT (IN) :: loplod 
+    INTEGER, INTENT (IN) :: nobd
     INTEGER, INTENT (IN) :: ne,jsp
     !     ..
     !     .. Array Arguments ..
@@ -56,7 +55,7 @@ CONTAINS
     !     .. Local Scalars ..
     INTEGER, PARAMETER :: lmaxb=3
     COMPLEX dtd,dtu,utd,utu
-    INTEGER lo 
+    INTEGER lo, mlotot, mlolotot, mlot_d, mlolot_d
     INTEGER i,ie,im,in,l1,l2,ll1,ll2,lm1,lm2,m1,m2,n,natom
     INTEGER natrun,is,isinv,j,irinv,it
     REAL   ,PARAMETER:: zero=0.0
@@ -73,11 +72,18 @@ CONTAINS
     !     ..
     !     ..
     !dimension%lmplmd = (dimension%lmd* (dimension%lmd+3))/2
+    mlotot = 0 ; mlolotot = 0
+    DO n = 1, atoms%ntype
+       mlotot = mlotot + atoms%nlo(n)
+       mlolotot = mlolotot + atoms%nlo(n)*(atoms%nlo(n)+1)/2
+    ENDDO
+    mlot_d = max(mlotot,1)
+    mlolot_d = max(mlolotot,1)
     ALLOCATE ( tlmplm%tdd(0:dimension%lmplmd,atoms%ntype,1),tlmplm%tuu(0:dimension%lmplmd,atoms%ntype,1),&
          tlmplm%tdu(0:dimension%lmplmd,atoms%ntype,1),tlmplm%tud(0:dimension%lmplmd,atoms%ntype,1),&
-         tlmplm%tuulo(0:dimension%lmd,-atoms%llod:atoms%llod,atoms%nlod,1),&
-         tlmplm%tdulo(0:dimension%lmd,-atoms%llod:atoms%llod,atoms%nlod,1),&
-         tlmplm%tuloulo(-atoms%llod:atoms%llod,-atoms%llod:atoms%llod,loplod,1),&
+         tlmplm%tuulo(0:dimension%lmd,-atoms%llod:atoms%llod,mlot_d,1),&
+         tlmplm%tdulo(0:dimension%lmd,-atoms%llod:atoms%llod,mlot_d,1),&
+         tlmplm%tuloulo(-atoms%llod:atoms%llod,-atoms%llod:atoms%llod,mlolot_d,1),&
          v_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb),&
          a21(3,atoms%natd),b4(3,atoms%natd),tlmplm%ind(0:dimension%lmd,0:dimension%lmd,atoms%ntype,1) )
     !
@@ -179,11 +185,13 @@ CONTAINS
           !
           !--->    add the local orbital and U contribution to a21
           !
-          CALL force_a21_lo(nobd,atoms, loplod,jsp,n,we,eig,ne,&
-               acof,bcof,ccof,aveccof,bveccof,cveccof, tlmplm,usdus, a21)
+          CALL force_a21_lo(nobd,atoms,jsp,n,we,eig,ne,&
+                            acof,bcof,ccof,aveccof,bveccof,&
+                            cveccof, tlmplm,usdus, a21)
 
-          CALL force_a21_U(nobd,atoms,lmaxb,n,jsp,we,ne, usdus,v_mmp,&
-               acof,bcof,ccof,aveccof,bveccof,cveccof, a21)
+          CALL force_a21_U(nobd,atoms,lmaxb,n,jsp,we,ne,&
+                           usdus,v_mmp,acof,bcof,ccof,&
+                           aveccof,bveccof,cveccof, a21)
 
 #ifdef CPP_APW
           ! -> B4 force

force/force_a21_lo.f90

 MODULE m_forcea21lo
 CONTAINS
-  SUBROUTINE force_a21_lo(nobd,atoms, loplod,isp,itype,we,eig,ne,&
-       acof,bcof,ccof,aveccof,bveccof,cveccof, tlmplm,usdus, a21)
+  SUBROUTINE force_a21_lo(nobd,atoms,isp,itype,we,eig,ne,&
+                          acof,bcof,ccof,aveccof,bveccof,&
+                          cveccof,tlmplm,usdus, a21)
     !
     !***********************************************************************
     ! This subroutine calculates the local orbital contribution to A21,
@@ -18,7 +19,7 @@ CONTAINS
     !     ..
     !     .. Scalar Arguments ..
     INTEGER, INTENT (IN) :: nobd     
-    INTEGER, INTENT (IN) :: loplod,itype,ne,isp
+    INTEGER, INTENT (IN) :: itype,ne,isp
     !     ..
     !     .. Array Arguments ..
     REAL,    INTENT (IN) :: we(nobd),eig(:)!(dimension%neigd)

force/force_w.f90

@@ -8,6 +8,7 @@
       USE m_geo
       USE m_relax
       USE m_types
+      USE m_xmlOutput
       IMPLICIT NONE
 
       TYPE(t_results),INTENT(IN)   :: results
@@ -26,6 +27,7 @@
 !     ..
 !     .. Local Arrays ..
       REAL forcetot(3,atoms%ntypd)
+      CHARACTER(LEN=20) :: attributes(7)
 !
 !     write spin-dependent forces
 !
@@ -52,30 +54,43 @@
       WRITE  (6,8005)
       WRITE (16,8005)
  8005 FORMAT (/,' ***** TOTAL FORCES ON ATOMS ***** ',/)
+      IF (input%l_f) CALL openXMLElement('totalForcesOnRepresentativeAtoms',(/'units'/),(/'Htr/bohr'/))
       nat1 = 1
       DO n = 1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
-!
-         DO i = 1,3
-            forcetot(i,n) = zero
-         END DO
-         DO jsp = 1,input%jspins
             DO i = 1,3
-               forcetot(i,n) = forcetot(i,n) + results%force(i,n,jsp)
+               forcetot(i,n) = zero
             END DO
-         END DO
-!
-         WRITE (6,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
-     &     (forcetot(i,n),i=1,3)
-         WRITE (16,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
-     &     (forcetot(i,n),i=1,3)
- 8010    FORMAT (' TOTAL FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,3x,'Y=',&
-     &          f7.3,3x,'Z=',f7.3,/,22x,' FX_TOT=',f9.6,&
-     &          ' FY_TOT=',f9.6,' FZ_TOT=',f9.6)
-!
-         ENDIF
+            DO jsp = 1,input%jspins
+               DO i = 1,3
+                  forcetot(i,n) = forcetot(i,n) + results%force(i,n,jsp)
+               END DO
+            END DO
+
+            WRITE (6,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
+     &        (forcetot(i,n),i=1,3)
+            WRITE (16,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
+     &        (forcetot(i,n),i=1,3)
+ 8010       FORMAT (' TOTAL FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,3x,'Y=',&
+     &              f7.3,3x,'Z=',f7.3,/,22x,' FX_TOT=',f9.6,&
+     &              ' FY_TOT=',f9.6,' FZ_TOT=',f9.6)
+
+            WRITE(attributes(1),'(i0)') n
+            WRITE(attributes(2),'(f12.6)') atoms%pos(1,nat1)
+            WRITE(attributes(3),'(f12.6)') atoms%pos(2,nat1)
+            WRITE(attributes(4),'(f12.6)') atoms%pos(3,nat1)
+            WRITE(attributes(5),'(f12.8)') forcetot(1,n)
+            WRITE(attributes(6),'(f12.8)') forcetot(2,n)
+            WRITE(attributes(7),'(f12.8)') forcetot(3,n)
+            IF (input%l_f) THEN
+               CALL writeXMLElementFormPoly('forceTotal',(/'atomType','x       ','y       ','z       ',&
+                                                           'F_x     ','F_y     ','F_z     '/),attributes,&
+                                            reshape((/8,1,1,1,3,3,3,6,12,12,12,12,12,12/),(/7,2/)))
+            END IF
+         END IF
          nat1 = nat1 + atoms%neq(n)
       END DO
+      IF (input%l_f) CALL closeXMLElement('totalForcesOnRepresentativeAtoms')
 
       sum=0.0
       DO n = 1,atoms%ntype