Adding MARTINI
[unres4.git] / source / cluster / io_clust.F90
index ed4ccc3..312b9f8 100644 (file)
       call int_from_cart1(.false.)
       do j=nnt+1,nct-1
         mnum=molnum(j)
-!        write (iout,*) "Check atom",j
+        write (iout,*) "Check atom",j,itype(j,mnum),vbld(j)
         if (mnum.ne.1) cycle
         if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
-        if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
+        if (itype(j-1,molnum(j-1)).eq.ntyp1_molec(molnum(j-1))) cycle
 
         if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
          if (j.gt.2) then
           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
-          write (iout,*) "Conformation",jjj,jj+1
+          write (iout,*) "Conformation",jjj,jj+1,itype(j-1,molnum(j-1)),itype(j,mnum)
           write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
        chalen
           write (iout,*) "The Cartesian geometry is:"
       use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
                  iscode,symetr,punch_dist,print_dist,nstart,nend,&
                  caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
-                 r_cut_ele,nclust,tor_mode,scelemode
+                 r_cut_ele,nclust,tor_mode,scelemode,r_cut_mart,r_cut_ang,&
+                 rlamb_mart
+      use geometry_data, only:bordliptop,bordlipbot,&
+          bufliptop,buflipbot,lipthick,lipbufthick
 !      implicit none
 !      include 'DIMENSIONS'
 !      include 'sizesclu.dat'
       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+      call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+      call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+      if (lipthick.gt.0.0d0) then
+       bordliptop=(boxzsize+lipthick)/2.0
+       bordlipbot=bordliptop-lipthick
+      if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+      write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+      buflipbot=bordlipbot+lipbufthick
+      bufliptop=bordliptop-lipbufthick
+      if ((lipbufthick*2.0d0).gt.lipthick) &
+       write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+      endif !lipthick.gt.0
+      write(iout,*) "bordliptop=",bordliptop
+      write(iout,*) "bordlipbot=",bordlipbot
+      write(iout,*) "bufliptop=",bufliptop
+      write(iout,*) "buflipbot=",buflipbot
+
       call readi(controlcard,'NCLUST',nclust,5)
 !      ions=index(controlcard,"IONS").gt.0
       call readi(controlcard,"SCELEMODE",scelemode,0)
       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
+      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)
       call readi(controlcard,'NEND',nend,0)
       call reada(controlcard,'ECUT',ecut,10.0d0)
       call reada(controlcard,'PROB',prob_limit,0.99d0)
       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
-
+      call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
+      call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
       call reada(weightcard,"D0CM",d0cm,3.78d0)
       call reada(weightcard,"AKCM",akcm,15.1d0)
       call reada(weightcard,"AKTH",akth,11.0d0)
       call reada(weightcard,"DTRISS",dtriss,1.001D0)
       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
+      call reada(weightcard,"LIPSCALE",lipscale,1.0D0)
+      call reada(weightcard,"WTL",wliptran,1.0D0)
 
       call reada(weightcard,"HT",Ht,0.0D0)
       if (dyn_ss) then
           weights(46)=wscbase
           weights(47)=wscpho
           weights(48)=wpeppho
-
+          weights(49)=wpeppho
+          weights(50)=wcatnucl
+          weights(56)=wcat_tran
+          weights(22)=wliptran
 
       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
@@ -1870,6 +1899,11 @@ write(iout,*)"po setup var"
       open (itube,file=tubename,status='old')
       call getenv('IONPAR',ionname)
       open (iion,file=ionname,status='old')
+      call getenv('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old')
+      call getenv('IONPAR_TRAN',iontranname)
+      open (iiontran,file=iontranname,status='old')
+
 
 #ifndef OLDSCP
 !
@@ -1904,7 +1938,7 @@ write(iout,*)"po setup var"
       character(len=50) :: tytul
       character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
                                                'G','H','I','J'/)
-      integer :: ica(nres),k,iti1
+      integer :: ica(nres),k,iti1,ki
       real(kind=8) :: etot,rmsd
       integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
       real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
@@ -1946,11 +1980,21 @@ write(iout,*)"po setup var"
           write (ipdb,'(a)') 'TER'
         else
         ires=ires+1
-        if (ires.eq.2) then
-         do j=1,3
-         cbeg(j)=c(j,i-1)
-         enddo
-        endif
+!         if (molnum(i).ge.1) then
+          if (i.le.3) then
+           ki=2
+          else
+          if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
+          ki=ki+1
+          endif
+         endif
+!        endif
+!          print *,"tu2",i,k
+          if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
+          if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
+        do j=1,3
+        cbeg(j)=c(j,ki)
+        enddo
         iatom=iatom+1
         ica(i)=iatom
         call to_box(c(1,i),c(2,i),c(3,i))