Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / unres / io_config.F90
index 44ec36c..4fc0905 100644 (file)
                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
                 res1,epsijlip,epspeptube,epssctube,sigmapeptube,      &
-                sigmasctube
-      integer :: ichir1,ichir2,ijunk
+                sigmasctube,krad2,ract
+      integer :: ichir1,ichir2,ijunk,irdiff
       character*3 string
-
+      character*80 temp1,mychar
 !      real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
 !el      allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 !el      allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
       msc(:,:)=0.0d0
       isc(:,:)=0.0d0
 
-      allocate(msc(ntyp+1,5)) !(ntyp+1)
-      allocate(isc(ntyp+1,5)) !(ntyp+1)
-      allocate(restok(ntyp+1,5)) !(ntyp+1)
+      allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
+      allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
+      allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1)
 
       read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
       do i=1,ntyp_molec(1)
 
 
 
-      if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
+      if (.not.allocated(ichargecat)) &
+      allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
+      ichargecat(:)=0
        if (oldion.eq.1) then
             do i=1,ntyp_molec(5)
              read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
         print *,liptranene(i)
        enddo
        close(iliptranpar)
+! water parmaters entalphy
+       allocate(awaterenta(0:400))
+       allocate(bwaterenta(0:400))
+       allocate(cwaterenta(0:400))
+       allocate(dwaterenta(0:400))
+       allocate(awaterentro(0:400))
+       allocate(bwaterentro(0:400))
+       allocate(cwaterentro(0:400))
+       allocate(dwaterentro(0:400))
+
+       read(iwaterwater,*) mychar
+       read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
+                           cwaterenta(0),dwaterenta(0)
+       do i=1,398
+       read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
+                           cwaterenta(i),dwaterenta(i)
+       irdiff=int((ract-2.06d0)*50.0d0)+1
+       if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
+       enddo
+! water parmaters entrophy
+       read(iwaterwater,*) mychar
+       read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
+                           cwaterentro(0),dwaterentro(0)
+       do i=1,398
+       read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
+                           cwaterentro(i),dwaterentro(i)
+       irdiff=int((ract-2.06d0)*50.0d0)+1
+       if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
+       enddo
+
 
 #ifdef CRYST_THETA
 !
       allocate(ddnewtor(3,2,-nloctyp:nloctyp))
       allocate(e0newtor(3,-nloctyp:nloctyp))
       allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
-
+      bnew1=0.0d0
+      bnew2=0.0d0
+      ccnew=0.0d0
+      ddnew=0.0d0
+      e0new=0.0d0
+      eenew=0.0d0
+      bnew1tor=0.0d0
+      bnew2tor=0.0d0
+      ccnewtor=0.0d0
+      ddnewtor=0.0d0
+      e0newtor=0.0d0
+      eenewtor=0.0d0
       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
       itype2loc(ntyp1)=nloctyp
        print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
        END DO
       END DO
+        do i=1,ntyp
+         read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
+        enddo
       do i=1,ntyp
        do j=1,i
         IF ((LaTeX).and.(i.gt.24)) then
       allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
       allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
 
-
+      write (iout,*) "ESCBASEPARM"
       do i=1,ntyp_molec(1)
-       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
-        write (*,*) "Im in ", i, " ", j
+       do j=1,ntyp_molec(2) ! without U then we will take T for U
+!        write (*,*) "Im in ", i, " ", j
         read(isidep_scbase,*) &
         eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
         chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
-         write(*,*) "eps",eps_scbase(i,j)
+!         write(*,*) "eps",eps_scbase(i,j)
         read(isidep_scbase,*) &
        (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
        chis_scbase(i,j,1),chis_scbase(i,j,2)
         read(isidep_scbase,*) &
        alphapol_scbase(i,j), &
        epsintab_scbase(i,j) 
+        if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
+        if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
+        if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
+        if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
+        if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
+        if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
+        if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
+        if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
+        write(iout,*) &
+        eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+        chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+         write(*,*) "eps",eps_scbase(i,j)
+        write(iout,*) &
+       (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+       chis_scbase(i,j,1),chis_scbase(i,j,2)
+        write(iout,*) &
+       dhead_scbasei(i,j), &
+       dhead_scbasej(i,j), &
+       rborn_scbasei(i,j),rborn_scbasej(i,j)
+        write(iout,*) &
+       (wdipdip_scbase(k,i,j),k=1,3), &
+       (wqdip_scbase(k,i,j),k=1,2)
+        write(iout,*) &
+       alphapol_scbase(i,j), &
+       epsintab_scbase(i,j)
+
        END DO
+        j=4
+        write(iout,*) &
+        eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+        chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+         write(*,*) "eps",eps_scbase(i,j)
+        write(iout,*) &
+       (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+       chis_scbase(i,j,1),chis_scbase(i,j,2)
+        write(iout,*) &
+       dhead_scbasei(i,j), &
+       dhead_scbasej(i,j), &
+       rborn_scbasei(i,j),rborn_scbasej(i,j)
+        write(iout,*) &
+       (wdipdip_scbase(k,i,j),k=1,3), &
+       (wqdip_scbase(k,i,j),k=1,2)
+        write(iout,*) &
+       alphapol_scbase(i,j), &
+       epsintab_scbase(i,j)
+
       END DO
       allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
       allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
       allocate(chis_pepbase(ntyp_molec(2),2))
       allocate(wdipdip_pepbase(3,ntyp_molec(2)))
 
+      write (iout,*) "EPEPBASEPARM"
 
-       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+       do j=1,ntyp_molec(2) ! without U then we will take T for U
         write (*,*) "Im in ", i, " ", j
         read(isidep_pepbase,*) &
         eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
         chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+        if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
+        if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
+        if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
+        if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
+        if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
+        if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
+        if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
+        if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
+
          write(*,*) "eps",eps_pepbase(j)
         read(isidep_pepbase,*) &
        (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
        chis_pepbase(j,1),chis_pepbase(j,2)
         read(isidep_pepbase,*) &
        (wdipdip_pepbase(k,j),k=1,3)
+        write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+        chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+        write(iout,*) &
+       (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+       chis_pepbase(j,1),chis_pepbase(j,2)
+        write(iout,*) &
+       (wdipdip_pepbase(k,j),k=1,3)
+
        END DO
+        j=4
+        write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+        chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+        write(iout,*) &
+       (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+       chis_pepbase(j,1),chis_pepbase(j,2)
+        write(iout,*) &
+       (wdipdip_pepbase(k,j),k=1,3)
+
       allocate(aa_pepbase(ntyp_molec(2)))
       allocate(bb_pepbase(ntyp_molec(2)))
 
         read(isidep_scpho,*) &
          epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
          alphi_scpho(j)
+        if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
+        if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
+        if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
+        if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
+        if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
+        if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
+        if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
+        if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
+
        
        END DO
       allocate(aa_scpho(ntyp_molec(1)))
 !      v1ss=0.0d0
 !      v2ss=0.0d0
 !      v3ss=0.0d0
+! MARTINI PARAMETER
+      allocate(ichargelipid(ntyp_molec(4)))
+      allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
+      allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
+      allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
+      allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
+      allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
+      kjtokcal=0.2390057361
+      krad=57.295779513
+      !HERE THE MASS of MARTINI
+      write(*,*) "before MARTINI PARAM"
+      do i=1,ntyp_molec(4)
+       msc(i,4)=0.0d0
+       mp(4)=72.0d0
+       isc(i,4)=0.d0
+      enddo
+      ip(4)=0.0
+      !relative dielectric constant = 15 for implicit screening
+      k_coulomb_lip=332.0d0/15.0d0
+      !kbond = 1250 kJ/(mol*nm*2)
+      kbondlip=1250.0d0*kjtokcal/100.0d0
+      krad2=krad**2.0
+      lip_angle_force=0.0d0
+      if (DRY_MARTINI.gt.0) then
+      lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
+      lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
+      lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
+      lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
+      lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
+      lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
+      lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
+      lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
+      else
+      lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
+      lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
+      lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
+      lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
+      lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
+      lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
+      lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
+      lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
+      endif
+      lip_angle_angle=0.0d0
+      lip_angle_angle(3,12,12)=120.0/krad
+      lip_angle_angle(3,12,18)=180.0/krad
+      lip_angle_angle(3,18,16)=180.0/krad
+      lip_angle_angle(12,18,16)=180.0/krad
+      lip_angle_angle(18,16,18)=120.0/krad
+      lip_angle_angle(16,18,18)=180.0/krad
+      lip_angle_angle(12,18,18)=180.0/krad
+      lip_angle_angle(18,18,18)=180.0/krad
+       read(ilipbond,*) temp1
+      do i=1,18
+        read(ilipbond,*) temp1, lip_bond(i,1), &
+        lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
+        lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
+        lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
+        lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
+        lip_bond(i,18)
+        do j=1,18
+          lip_bond(i,j)=lip_bond(i,j)*10
+        enddo
+      enddo
+
+       read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
+       read(ilipnonbond,*) temp1
+      do i=1,18
+        read(ilipnonbond,*) temp1, lip_eps(i,1), &
+        lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
+        lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
+        lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
+        lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
+        lip_eps(i,18)
+!        write(*,*) i, lip_eps(i,18)
+        do j=1,18
+          lip_eps(i,j)=lip_eps(i,j)*kjtokcal
+        enddo
+      enddo
+       read(ilipnonbond,*) temp1
+      do i=1,18
+        read(ilipnonbond,*) temp1,lip_sig(i,1), &
+        lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
+        lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
+        lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
+        lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
+        lip_sig(i,18)
+        do j=1,18
+          lip_sig(i,j)=lip_sig(i,j)*10.0
+        enddo
+      enddo
+      write(*,*) "after MARTINI PARAM"
 
 ! Ions by Aga
 
       if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
 
 
-            if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
+            if (.not.allocated(ichargecat))&
+       allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
+      write(*,*) "before ions",oldion
+            ichargecat(:)=0
+
 ! 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
             else
              read(iion,*) ijunk
             endif
-
-            do i=1,ntyp_molec(5)
+            print *,ntyp_molec(5)
+            do i=-ntyp_molec(5),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
            ! mp(5)=0.2
-           ! pstok(5)=1.0
+             pstok(5)=3.0
 !DIR$ NOUNROLL 
-      do j=1,ntyp_molec(5)
+      do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
        do i=1,ntyp
 !       do j=1,ntyp_molec(5)
 !        write (*,*) "Im in ALAB", i, " ", j
        enddo
 
        do i=1,ntyp
-       do j=1,ntyp_molec(5)
+       do j=1,ntyp_molec(5)-1
       if (i.eq.10) then
       write (iout,*) 'i= ', i, ' j= ', j
       write (iout,*) 'epsi0= ', epscat(i,j)
        enddo
 
       endif
-
-      
+! read number of Zn2+
+! here two denotes the Zn2+ and Cu2+
+      write(iout,*) "before TRANPARM"
+      allocate(aomicattr(0:3,2))
+      allocate(athetacattran(0:6,5,2))
+      allocate(agamacattran(3,5,2))
+      allocate(acatshiftdsc(5,2))
+      allocate(bcatshiftdsc(5,2))
+      allocate(demorsecat(5,2))
+      allocate(alphamorsecat(5,2))
+      allocate(x0catleft(5,2))
+      allocate(x0catright(5,2))
+      allocate(x0cattrans(5,2))
+      allocate(ntrantyp(2))
+      do i=1,1 ! currently only Zn2+
+
+      read(iiontran,*) ntrantyp(i)
+!ntrantyp=4
+!| ao0 ao1 ao2 ao3
+!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
+!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
+!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+      read(iiontran,*) (aomicattr(j,i),j=0,3)
+      do j=1,ntrantyp(i)
+       read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
+       (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+       demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+       x0cattrans(j,i)
+      enddo 
+      enddo
       if(me.eq.king) then
       write (iout,'(/a)') "Disulfide bridge parameters:"
       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
   119 write (iout,*) "Error reading SCCOR parameters"
       go to 999
   121 write (iout,*) "Error in Czybyshev parameters"
+      go to 999
+  123 write(iout,*) "Error in transition metal parameters"
   999 continue
 #ifdef MPI
       call MPI_Finalize(Ierror)
           itype(ires_old,molecule)=ntyp1_molec(molecule)
           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
           nres_molec(molecule)=nres_molec(molecule)+2
+!          if (molecule.eq.4) ires=ires+2
           ibeg=2
 !          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
 !          iii=0
           endif
           iii=0
+        else if (card(:3).eq.'BRA') then
+         molecule=4
+          ires=ires+1
+          ires_old=ires+1
+          itype(ires,molecule)=ntyp1_molec(molecule)-1
+          nres_molec(molecule)=nres_molec(molecule)+1
+        
         endif
 ! Read free energy
         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
           read (card(12:16),*) atom
 !          write (iout,*) "! ",atom," !",ires
 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+         if (card(14:16).eq.'LIP') then
+! reading lipid
+          if (ibeg.eq.1) then
+          molecule=4
+          ires=ires+1
+                nres_molec(molecule)=nres_molec(molecule)+1
+                   itype(ires,molecule)=ntyp1_molec(molecule)
+          ibeg=0
+          endif
+         if (ibeg.eq.2) then
+         ibeg=0
+         ires=ires+2 
+         endif
+
+          molecule=4 
+                nres_molec(molecule)=nres_molec(molecule)+1
+          read (card(18:20),'(a3)') res
+          ires=ires+1
+          read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+          if (UNRES_PDB) then!
+              itype(ires,molecule)=rescode(ires,res,0,molecule)
+          else
+             itype(ires,molecule)=rescode_lip(res,ires)
+          endif
+         else 
           read (card(23:26),*) ires
           read (card(18:20),'(a3)') res
 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
               endif
           endif
+          endif !LIP
 !         print *,"IONS",ions,card(1:6)
         else if ((ions).and.(card(1:6).eq.'HETATM')) then
        if (firstion.eq.0) then 
        endif ! unres_pdb
        endif !firstion
           read (card(12:16),*) atom
-!          print *,"HETATOM", atom
+          print *,"HETATOM", atom(1:2)
           read (card(18:20),'(a3)') res
+          if (atom(3:3).eq.'H') cycle
           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
-          .or.(atom(1:2).eq.'K ')) &
-          then
+          .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
+          (atom(1:2).eq.'O '))  then
            ires=ires+1
+           print *,"I have water"
            if (molecule.ne.5) molecprev=molecule
            molecule=5
            nres_molec(molecule)=nres_molec(molecule)+1
            print *,"HERE",nres_molec(molecule)
-           res=res(2:3)//' '
+           if (res.ne.'WAT')  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! NA
 !      print *,"molecule",molecule
       if ((itype(nres,1).ne.10)) then
         nres=nres+1
+        nsup=nsup+1
           if (molecule.eq.5) molecule=molecprev
         itype(nres,molecule)=ntyp1_molec(molecule)
         nres_molec(molecule)=nres_molec(molecule)+1
           do j=1,3
           c_temporary(j,seqalingbegin)=c(j,i)
           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
-
           enddo
           itype_temporary(seqalingbegin,k)=itype(i,k)
           print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
        do i=1,2*nres
         do j=1,3
         c(j,i)=c_temporary(j,i)
+        if (i.gt.nres) then
+        if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
+        endif
         enddo
        enddo
        do k=1,5
        molnum(1)=5
        itype(1,5)=itype(2,5)
        itype(1,1)=0
-       do i=2,nres
+       do j=1,3
+        c(j,1)=c(j,2)
+       enddo
+       do i=2,nres-1
          itype(i,5)=itype(i+1,5)
+         do j=1,3
+          c(j,i)=c(j,i+1)
+         enddo
        enddo
        itype(nres,5)=0
        nres=nres-1
       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
       protein=index(controlcard,"PROTEIN").gt.0
       ions=index(controlcard,"IONS").gt.0
+      fodson=index(controlcard,"FODSON").gt.0
       call readi(controlcard,'OLDION',oldion,1)
       nucleic=index(controlcard,"NUCLEIC").gt.0
       write (iout,*) "with_theta_constr ",with_theta_constr
        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+       call reada(controlcard,"VNANO",velnanoconst,0.0d0)
        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
       write(iout,*) "R_CUT_ELE=",r_cut_ele
-! Lipidic parameters
+      call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
+      call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
+      call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
+
+! Lipidec parameters
       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
       if (lipthick.gt.0.0d0) then
 !
 ! Read MD settings
 !
-      use control_data, only: r_cut,rlamb,out1file
+      use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
       use energy_data
       use geometry_data, only: pi
       use MPI_data
       call readi(controlcard,"NSTEP",n_timestep,1000000)
       call readi(controlcard,"NTWE",ntwe,100)
       call readi(controlcard,"NTWX",ntwx,1000)
+      call readi(controlcard,"NFOD",nfodstep,100)
       call reada(controlcard,"DT",d_time,1.0d-1)
       call reada(controlcard,"DVMAX",dvmax,2.0d1)
       call reada(controlcard,"DAMAX",damax,1.0d1)
         enddo
 
         if(me.eq.king.or..not.out1file)then
+         do j=1,5
          write (iout,'(/2a/)') &
          "Radii of site types and friction coefficients and std's of",&
          " stochastic forces of fully exposed sites"
-         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
+         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
+        
          do i=1,ntyp
-          write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
-           gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
+          write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
+           gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
          enddo
+        enddo
         endif
       else if (tbf) then
         if(me.eq.king.or..not.out1file)then
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('IONPAR',ionname)
       open (iion,file=ionname,status='old',action='read')
-
+      call getenv_loc('IONPAR_TRAN',iontranname)
+      open (iiontran,file=iontranname,status='old',action='read')
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
 #elif (defined G77)
       open (iion,file=ionname,status='old')
       call getenv_loc('IONPAR_NUCL',ionnuclname)
       open (iionnucl,file=ionnuclname,status='old')
+      call getenv_loc('IONPAR_TRAN',iontranname)
+      open (iiontran,file=iontranname,status='old')
+      call getenv_loc('WATWAT',iwaterwatername)
+      open (iwaterwater,file=iwaterwatername,status='old')
+      call getenv_loc('WATPROT',iwaterscname)
+      open (iwatersc,file=iwaterscname,status='old')
+
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
 
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
+      call getenv_loc('LIPBOND',lipbondname)
+      open (ilipbond,file=lipbondname,status='old',action='read')
+      call getenv_loc('LIPNONBOND',lipnonbondname)
+      open (ilipnonbond,file=lipnonbondname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('IONPAR',ionname)
       open (iion,file=ionname,status='old',action='read')
       call getenv_loc('IONPAR_NUCL',ionnuclname)
       open (iionnucl,file=ionnuclname,status='old',action='read')
+      call getenv_loc('IONPAR_TRAN',iontranname)
+      open (iiontran,file=iontranname,status='old',action='read')
+      call getenv_loc('WATWAT',iwaterwatername)
+      open (iwaterwater,file=iwaterwatername,status='old',action='read')
+      call getenv_loc('WATPROT',iwaterscname)
+      open (iwatersc,file=iwaterscname,status='old',action='read')
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)