working gradient for cations
[unres4.git] / source / unres / io_config.F90
index a31dead..2906fb5 100644 (file)
           enddo
         enddo
       endif
+       if (oldion.eq.1) then
             do i=1,ntyp_molec(5)
              read(iion,*) msc(i,5),restok(i,5)
              print *,msc(i,5),restok(i,5)
             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
             enddo
             print *, catprm
+         endif
 !            read (iion,*) (vcatprm(k),k=1,ncatprotpram)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
         itype2loc(-i)=-itype2loc(i)
       enddo
 #else
+      allocate(iloctyp(-nloctyp:nloctyp))
+      allocate(itype2loc(-ntyp1:ntyp1))
       iloctyp(0)=10
       iloctyp(1)=9
       iloctyp(2)=20
       allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
       allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
       allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
-
+      allocate(b(13,-nloctyp-1:nloctyp+1))
       if (lprint) &
        write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
       do i=0,nloctyp-1
         b(4,-i)=-b(4,i)
         b(5,-i)=-b(5,i)
         endif
-c        B1(1,i)  = b(3)
-c        B1(2,i)  = b(5)
-c        B1(1,-i) = b(3)
-c        B1(2,-i) = -b(5)
-c        b1(1,i)=0.0d0
-c        b1(2,i)=0.0d0
-c        B1tilde(1,i) = b(3)
-c        B1tilde(2,i) =-b(5)
-c        B1tilde(1,-i) =-b(3)
-c        B1tilde(2,-i) =b(5)
-c        b1tilde(1,i)=0.0d0
-c        b1tilde(2,i)=0.0d0
-c        B2(1,i)  = b(2)
-c        B2(2,i)  = b(4)
-c        B2(1,-i)  =b(2)
-c        B2(2,-i)  =-b(4)
-cc        B1tilde(1,i) = b(3,i)
-cc        B1tilde(2,i) =-b(5,i)
-C        B1tilde(1,-i) =-b(3,i)
-C        B1tilde(2,-i) =b(5,i)
-cc        b1tilde(1,i)=0.0d0
-cc        b1tilde(2,i)=0.0d0
-cc        B2(1,i)  = b(2,i)
-cc        B2(2,i)  = b(4,i)
-C        B2(1,-i)  =b(2,i)
-C        B2(2,-i)  =-b(4,i)
-
-c        b2(1,i)=0.0d0
-c        b2(2,i)=0.0d0
+!c        B1(1,i)  = b(3)
+!c        B1(2,i)  = b(5)
+!c        B1(1,-i) = b(3)
+!c        B1(2,-i) = -b(5)
+!c        b1(1,i)=0.0d0
+!c        b1(2,i)=0.0d0
+!c        B1tilde(1,i) = b(3)
+!c!        B1tilde(2,i) =-b(5)
+!c!        B1tilde(1,-i) =-b(3)
+!c!        B1tilde(2,-i) =b(5)
+!c!        b1tilde(1,i)=0.0d0
+!c        b1tilde(2,i)=0.0d0
+!c        B2(1,i)  = b(2)
+!c        B2(2,i)  = b(4)
+!c        B2(1,-i)  =b(2)
+!c        B2(2,-i)  =-b(4)
+!cc        B1tilde(1,i) = b(3,i)
+!cc        B1tilde(2,i) =-b(5,i)
+!c        B1tilde(1,-i) =-b(3,i)
+!c        B1tilde(2,-i) =b(5,i)
+!cc        b1tilde(1,i)=0.0d0
+!cc        b1tilde(2,i)=0.0d0
+!cc        B2(1,i)  = b(2,i)
+!cc        B2(2,i)  = b(4,i)
+!c        B2(1,-i)  =b(2,i)
+!c        B2(2,-i)  =-b(4,i)
+
+!c        b2(1,i)=0.0d0
+!c        b2(2,i)=0.0d0
         CCold(1,1,i)= b(7,i)
         CCold(2,2,i)=-b(7,i)
         CCold(2,1,i)= b(9,i)
@@ -1878,27 +1882,27 @@ c        b2(2,i)=0.0d0
         CCold(2,2,-i)=-b(7,i)
         CCold(2,1,-i)=-b(9,i)
         CCold(1,2,-i)=-b(9,i)
-c        CC(1,1,i)=0.0d0
-c        CC(2,2,i)=0.0d0
-c        CC(2,1,i)=0.0d0
-c        CC(1,2,i)=0.0d0
-c        Ctilde(1,1,i)= CCold(1,1,i)
-c        Ctilde(1,2,i)= CCold(1,2,i)
-c        Ctilde(2,1,i)=-CCold(2,1,i)
-c        Ctilde(2,2,i)=-CCold(2,2,i)
-c        CC(1,1,i)=0.0d0
-c        CC(2,2,i)=0.0d0
-c        CC(2,1,i)=0.0d0
-c        CC(1,2,i)=0.0d0
-c        Ctilde(1,1,i)= CCold(1,1,i)
-c        Ctilde(1,2,i)= CCold(1,2,i)
-c        Ctilde(2,1,i)=-CCold(2,1,i)
-c        Ctilde(2,2,i)=-CCold(2,2,i)
-
-c        Ctilde(1,1,i)=0.0d0
-c        Ctilde(1,2,i)=0.0d0
-c        Ctilde(2,1,i)=0.0d0
-c        Ctilde(2,2,i)=0.0d0
+!c        CC(1,1,i)=0.0d0
+!c        CC(2,2,i)=0.0d0
+!c        CC(2,1,i)=0.0d0
+!c        CC(1,2,i)=0.0d0
+!c        Ctilde(1,1,i)= CCold(1,1,i)
+!c        Ctilde(1,2,i)= CCold(1,2,i)
+!c        Ctilde(2,1,i)=-CCold(2,1,i)
+!c        Ctilde(2,2,i)=-CCold(2,2,i)
+!c        CC(1,1,i)=0.0d0
+!c        CC(2,2,i)=0.0d0
+!c        CC(2,1,i)=0.0d0
+!c        CC(1,2,i)=0.0d0
+!c        Ctilde(1,1,i)= CCold(1,1,i)
+!c        Ctilde(1,2,i)= CCold(1,2,i)
+!c        Ctilde(2,1,i)=-CCold(2,1,i)
+!c        Ctilde(2,2,i)=-CCold(2,2,i)
+
+!c        Ctilde(1,1,i)=0.0d0
+!c        Ctilde(1,2,i)=0.0d0
+!c        Ctilde(2,1,i)=0.0d0
+!c        Ctilde(2,2,i)=0.0d0
         DDold(1,1,i)= b(6,i)
         DDold(2,2,i)=-b(6,i)
         DDold(2,1,i)= b(8,i)
@@ -1907,19 +1911,19 @@ c        Ctilde(2,2,i)=0.0d0
         DDold(2,2,-i)=-b(6,i)
         DDold(2,1,-i)=-b(8,i)
         DDold(1,2,-i)=-b(8,i)
-c        DD(1,1,i)=0.0d0
-c        DD(2,2,i)=0.0d0
-c        DD(2,1,i)=0.0d0
-c        DD(1,2,i)=0.0d0
-c        Dtilde(1,1,i)= DD(1,1,i)
-c        Dtilde(1,2,i)= DD(1,2,i)
-c        Dtilde(2,1,i)=-DD(2,1,i)
-c        Dtilde(2,2,i)=-DD(2,2,i)
-
-c        Dtilde(1,1,i)=0.0d0
-c        Dtilde(1,2,i)=0.0d0
-c        Dtilde(2,1,i)=0.0d0
-c        Dtilde(2,2,i)=0.0d0
+!c        DD(1,1,i)=0.0d0
+!c        DD(2,2,i)=0.0d0
+!c        DD(2,1,i)=0.0d0
+!c        DD(1,2,i)=0.0d0
+!c        Dtilde(1,1,i)= DD(1,1,i)
+!c        Dtilde(1,2,i)= DD(1,2,i)
+!c        Dtilde(2,1,i)=-DD(2,1,i)
+!c        Dtilde(2,2,i)=-DD(2,2,i)
+
+!c        Dtilde(1,1,i)=0.0d0
+!c        Dtilde(1,2,i)=0.0d0
+!c        Dtilde(2,1,i)=0.0d0
+!c        Dtilde(2,2,i)=0.0d0
         EEold(1,1,i)= b(10,i)+b(11,i)
         EEold(2,2,i)=-b(10,i)+b(11,i)
         EEold(2,1,i)= b(12,i)-b(13,i)
@@ -1930,11 +1934,11 @@ c        Dtilde(2,2,i)=0.0d0
         EEold(1,2,-i)=-b(12,i)-b(13,i)
         write(iout,*) "TU DOCHODZE"
         print *,"JESTEM"
-c        ee(1,1,i)=1.0d0
-c        ee(2,2,i)=1.0d0
-c        ee(2,1,i)=0.0d0
-c        ee(1,2,i)=0.0d0
-c        ee(2,1,i)=ee(1,2,i)
+!c        ee(1,1,i)=1.0d0
+!c        ee(2,2,i)=1.0d0
+!c        ee(2,1,i)=0.0d0
+!c        ee(1,2,i)=0.0d0
+!c        ee(2,1,i)=ee(1,2,i)
       enddo
       if (lprint) then
       write (iout,*)
@@ -2185,6 +2189,12 @@ c        ee(2,1,i)=ee(1,2,i)
       enddo
       enddo
       endif
+#ifndef NEWCORR
+      do i=1,ntyp1
+        itype2loc(i)=itortyp(i)
+      enddo
+#endif
+
       ELSE IF (TOR_MODE.eq.1) THEN
 
 !C read valence-torsional parameters
@@ -2644,25 +2654,47 @@ c        ee(2,1,i)=ee(1,2,i)
        allocate(wstate(4,ntyp,ntyp))
        allocate(dhead(2,2,ntyp,ntyp))
        allocate(nstate(ntyp,ntyp))
+       allocate(debaykap(ntyp,ntyp))
+
       if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
       if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+
       do i=1,ntyp
        do j=1,i
 !        write (*,*) "Im in ALAB", i, " ", j
         read(isidep,*) &
-       eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
-       (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
-       chis(i,j),chis(j,i), &
-       nstate(i,j),(wstate(k,i,j),k=1,4), &
-       dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
-       dtail(1,i,j),dtail(2,i,j), &
-       epshead(i,j),sig0head(i,j), &
-       rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
-       alphapol(i,j),alphapol(j,i), &
-       (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
+       eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
+       (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &           !6 w tej linii
+       chis(i,j),chis(j,i), &                                         !2 w tej linii
+       nstate(i,j),(wstate(k,i,j),k=1,4), &                           !5 w tej lini - 1 integer pierwszy
+       dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&  ! 4 w tej linii
+       dtail(1,i,j),dtail(2,i,j), &                                   ! 2 w tej lini
+       epshead(i,j),sig0head(i,j), &                                  ! 2 w tej linii
+       rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &       ! 5 w tej linii
+       alphapol(i,j),alphapol(j,i), &                                 ! 2 w tej linii
+       (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
+        IF ((LaTeX).and.(i.gt.24)) then
+        write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
+       eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
+       (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &           !6 w tej linii
+       chis(i,j),chis(j,i)                                            !2 w tej linii
+        endif
 !       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
        END DO
       END DO
+      do i=1,ntyp
+       do j=1,i
+        IF ((LaTeX).and.(i.gt.24)) then
+        write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
+       dhead(1,1,i,j),dhead(2,1,i,j),&  ! 2 w tej linii
+       dtail(1,i,j),dtail(2,i,j), &                                   ! 2 w tej lini
+       epshead(i,j),sig0head(i,j), &                                  ! 2 w tej linii
+       rborn(i,j),rborn(j,i), &       ! 3 w tej linii
+       alphapol(i,j),alphapol(j,i), &                                 ! 2 w tej linii
+       (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
+        endif
+       END DO
+      END DO
       DO i = 1, ntyp
        DO j = i+1, ntyp
         eps(i,j) = eps(j,i)
@@ -2695,10 +2727,13 @@ c        ee(2,1,i)=ee(1,2,i)
 
         wquad(i,j) = wquad(j,i)
         epsintab(i,j) = epsintab(j,i)
+        debaykap(i,j)=debaykap(j,i)
 !        if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
        END DO
       END DO
       endif
+
+
 !      goto 50
 !--------------------- GBV potential -----------------------------------
        case(5)
@@ -2790,7 +2825,7 @@ c        ee(2,1,i)=ee(1,2,i)
          sigeps=dsign(1.0D0,epsij)
          epsij=dabs(epsij)
          aa_aq(i,j)=epsij*rrij*rrij
-          print *,"ADASKO",epsij,rrij,expon
+!          print *,"ADASKO",epsij,rrij,expon
          bb_aq(i,j)=-sigeps*epsij*rrij
          aa_aq(j,i)=aa_aq(i,j)
          bb_aq(j,i)=bb_aq(i,j)
@@ -3199,6 +3234,72 @@ c        ee(2,1,i)=ee(1,2,i)
 !      v1ss=0.0d0
 !      v2ss=0.0d0
 !      v3ss=0.0d0
+
+! Ions by Aga
+
+       allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+       allocate(sigiso1cat(ntyp,ntyp),rborncat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+       allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+       allocate(chiscat(ntyp,ntyp),wquadcat(ntyp,ntyp),chippcat(ntyp,ntyp))
+       allocate(epsintabcat(ntyp,ntyp))
+       allocate(dtailcat(2,ntyp,ntyp))
+       allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
+       allocate(wqdipcat(2,ntyp,ntyp))
+       allocate(wstatecat(4,ntyp,ntyp))
+       allocate(dheadcat(2,2,ntyp,ntyp))
+       allocate(nstatecat(ntyp,ntyp))
+       allocate(debaykapcat(ntyp,ntyp))
+
+      if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
+      if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+      if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+      allocate (ichargecat(ntyp_molec(5)))
+! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
+       if (oldion.eq.0) then
+            if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
+            allocate(icharge(1:ntyp1))
+            read(iion,*) (icharge(i),i=1,ntyp)
+            else
+             read(iion,*) ijunk
+            endif
+
+            do i=1,ntyp_molec(5)
+             read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
+             print *,msc(i,5),restok(i,5)
+            enddo
+            ip(5)=0.2
+
+      do i=1,ntyp
+       do j=1,ntyp_molec(5)
+!        write (*,*) "Im in ALAB", i, " ", j
+        read(iion,*) &
+       epscat(i,j),sigmacat(i,j),chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
+       (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+       chiscat(i,j),chiscat(j,i), &
+       nstatecat(i,j),(wstatecat(k,i,j),k=1,4), &                           !5 w tej lini - 1 integer pierwszy
+       dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+       dtailcat(1,i,j),dtailcat(2,i,j), &
+       epsheadcat(i,j),sig0headcat(i,j), &
+!wdipcat = w1 , w2
+       rborncat(i,j),rborncat(j,i),(wqdipcat(k,i,j),k=1,2), &
+       alphapolcat(i,j),alphapolcat(j,i), &
+       (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
+!       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
+       END DO
+      END DO
+      allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+      do i=1,ntyp
+        do j=1,ntyp_molec(5)
+          epsij=epscat(i,j)
+          rrij=sigmacat(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_aq_cat(i,j)=epsij*rrij*rrij
+          bb_aq_cat(i,j)=-sigeps*epsij*rrij
+         enddo
+       enddo 
+      endif
+
       
       if(me.eq.king) then
       write (iout,'(/a)') "Disulfide bridge parameters:"
@@ -3407,6 +3508,8 @@ c        ee(2,1,i)=ee(1,2,i)
         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 ! Fish out the ATOM cards.
 !        write(iout,*) 'card',card(1:20)
+!        print *,"ATU ",card(1:6), CARD(3:6)
+!        print *,card
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
 !          write (iout,*) "! ",atom," !",ires
@@ -3428,9 +3531,9 @@ c        ee(2,1,i)=ee(1,2,i)
                 enddo
               else
                 call sccenter(ires_old,iii,sccor)
-              endif
+              endif !unres_pdb
               iii=0
-            endif
+            endif !ind_pdb
 ! Start new residue.
             if (res.eq.'Cl-' .or. res.eq.'Na+') then
               ires=ires_old
@@ -3442,7 +3545,7 @@ c        ee(2,1,i)=ee(1,2,i)
                 ishift=ishift-1
                 itype(1,1)=ntyp1
                 nres_molec(molecule)=nres_molec(molecule)+1
-              endif
+              endif ! Gly
               ires=ires-ishift+ishift1
               ires_old=ires
 !              write (iout,*) "ishift",ishift," ires",ires,&
@@ -3461,7 +3564,7 @@ c        ee(2,1,i)=ee(1,2,i)
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
               ires=ires-ishift+ishift1
               ires_old=ires
-            endif 
+            endif ! Na Cl
 !            print *,'atom',ires,atom
             if (res.eq.'ACE' .or. res.eq.'NHE') then
               itype(ires,1)=10
@@ -3485,11 +3588,11 @@ c        ee(2,1,i)=ee(1,2,i)
 !              print *,"ires=",ires,istype(ires)
 !            endif
 
-            endif
-            endif
+            endif ! atom.eq.CA
+            endif !ACE
           else
             ires=ires-ishift+ishift1
-          endif
+          endif !ires_old
 !          write (iout,*) "ires_old",ires_old," ires",ires
           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 !            ishift1=ishift1+1
@@ -3574,6 +3677,7 @@ c        ee(2,1,i)=ee(1,2,i)
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
               endif
           endif
+!         print *,"IONS",ions,card(1:6)
         else if ((ions).and.(card(1:6).eq.'HETATM')) then
        if (firstion.eq.0) then 
        firstion=1
@@ -3583,10 +3687,10 @@ c        ee(2,1,i)=ee(1,2,i)
          enddo
        else
           call sccenter(ires,iii,sccor)
-       endif
-       endif
+       endif ! unres_pdb
+       endif !firstion
           read (card(12:16),*) atom
-          print *,"HETATOM", atom
+!          print *,"HETATOM", atom
           read (card(18:20),'(a3)') res
           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
@@ -3600,7 +3704,7 @@ c        ee(2,1,i)=ee(1,2,i)
            res=res(2:3)//' '
            itype(ires,molecule)=rescode(ires,res,0,molecule)
            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-          endif
+          endif! NA
         endif !atom
       enddo
    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
@@ -3674,7 +3778,8 @@ c        ee(2,1,i)=ee(1,2,i)
 !              e2(3)=0.0d0
 !            endif
             do j=1,3
-              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+!              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+             c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
             enddo
 !          else !unres_pdb
            do j=1,3
@@ -3789,7 +3894,8 @@ c        ee(2,1,i)=ee(1,2,i)
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-1.9d0*e2(j)
+!            c(j,1)=c(j,2)-1.9d0*e2(j)
+             c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
           enddo
         else
         do j=1,3
@@ -4206,6 +4312,7 @@ c        ee(2,1,i)=ee(1,2,i)
       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
       protein=index(controlcard,"PROTEIN").gt.0
       ions=index(controlcard,"IONS").gt.0
+      call readi(controlcard,'OLDION',oldion,1)
       nucleic=index(controlcard,"NUCLEIC").gt.0
       write (iout,*) "with_theta_constr ",with_theta_constr
       AFMlog=(index(controlcard,'AFM'))
@@ -4223,6 +4330,8 @@ c        ee(2,1,i)=ee(1,2,i)
 ! elemode = 0 is orignal UNRES electrostatics
 ! elemode = 1 is "Momo" potentials in progress
 ! elemode = 2 is in development EVALD
+
+
       write (iout,*) TUBEmode,"TUBEMODE"
       if (TUBEmode.gt.0) then
        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)