X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fio_clust.F90;h=312b9f843ee17d556336adccecdbea95a7ff51a9;hb=0f5304c0d44e3ec4d6538de7d274c4b1a1930049;hp=6caef1c5a89a64afcf3e321aa3d60d6fdc0aa2dd;hpb=a280f55746ca567cfd02f1184da61a28fbe9ebc7;p=unres4.git diff --git a/source/cluster/io_clust.F90 b/source/cluster/io_clust.F90 index 6caef1c..312b9f8 100644 --- a/source/cluster/io_clust.F90 +++ b/source/cluster/io_clust.F90 @@ -559,7 +559,9 @@ allocate(nss_all(maxstr_proc)) !(maxstr_proc) allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc) allocate(iscore(maxconf)) !(maxconf) - + nss_all=0 + ihpb_all=0 + jhpb_all=0 #ifdef CHUJ ICON=1 @@ -741,8 +743,8 @@ ! call flush(iout) if (iret.eq.0) goto 101 call xdrfint(ixdrf, nss, iret) -! write (iout,*) "iret",iret -! write (iout,*) "nss",nss + write (iout,*) "iret",iret +! write (iout,*) "nss",nss,i,"TUTU" call flush(iout) if (iret.eq.0) goto 101 do k=1,nss @@ -750,10 +752,13 @@ if (iret.eq.0) goto 101 call xdrfint(ixdrf, jhpb(k), iret) if (iret.eq.0) goto 101 +! write(iout,*), ihpb(k),jhpb(k),"TUTU" enddo call xdrffloat(ixdrf,reini,iret) if (iret.eq.0) goto 101 +! write(iout,*) "TUTU", reini call xdrffloat(ixdrf,refree,iret) +! write(iout,*) "TUTU", refree if (iret.eq.0) goto 101 call xdrffloat(ixdrf,rmsdev,iret) if (iret.eq.0) goto 101 @@ -889,15 +894,15 @@ 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:" @@ -1134,6 +1139,8 @@ !------------------------------------------------------------------------------ subroutine daread_ccoords(istart_conf,iend_conf) + use energy_data, only: dyn_ss + ! implicit none ! include "DIMENSIONS" ! include "sizesclu.dat" @@ -1172,10 +1179,17 @@ write (iout,*) "Reading binary file, record",iii," ii",ii call flush(iout) #endif + if (dyn_ss) then + read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), & + ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),& + entfac(ii),rmstb(ii) + else read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),& - ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),& - nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),& - entfac(ii),rmstb(ii) + ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),& + nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),& + entfac(ii),rmstb(ii) + endif + #ifdef DEBUG write (iout,*) ii,iii,ij,entfac(ii) write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres) @@ -1195,6 +1209,8 @@ ! implicit none ! include "DIMENSIONS" ! include "sizesclu.dat" + use energy_data, only: dyn_ss + #ifdef MPI use MPI_data include "mpif.h" @@ -1231,10 +1247,17 @@ write (iout,*) "Writing binary file, record",iii," ii",ii call flush(iout) #endif + if (dyn_ss) then write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),& - ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),& - nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),& - entfac(ii),rmstb(ii) + ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),& + entfac(ii),rmstb(ii) + else + write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), & + ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),& + nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),& + entfac(ii),rmstb(ii) + endif + #ifdef DEBUG write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres) write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,& @@ -1258,7 +1281,10 @@ 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' @@ -1316,6 +1342,23 @@ 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) @@ -1323,9 +1366,12 @@ call readi(controlcard,'TORMODE',tor_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "torsional and valence angle mode",tor_mode - call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.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 + 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) @@ -1390,6 +1436,7 @@ character(len=800) :: weightcard ! integer rescode real(kind=8) :: x(6*nres) !(maxvar) + real(kind=8) :: ssscale integer :: itype_pdb(nres) !(maxres) ! logical seq_comp integer :: i,j,kkk,mnum @@ -1469,7 +1516,38 @@ 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,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + call reada(weightcard,"V2SS",v2ss,7.61d0) + call reada(weightcard,"V3SS",v3ss,13.7d0) + call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"ATRISS",atriss,0.301D0) + call reada(weightcard,"BTRISS",btriss,0.021D0) + call reada(weightcard,"CTRISS",ctriss,1.001D0) + 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 + ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale + Ht=(Ht/wsc-0.25*eps(1,1))*ssscale + akcm=akcm/wsc*ssscale + akth=akth/wsc*ssscale + akct=akct/wsc*ssscale + v1ss=v1ss/wsc*ssscale + v2ss=v2ss/wsc*ssscale + v3ss=v3ss/wsc*ssscale + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc @@ -1510,7 +1588,10 @@ 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,& @@ -1672,6 +1753,23 @@ write(iout,*)"po ini_int_tab" write(iout,*)"przed setup var" call setup_var write(iout,*)"po setup var" + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + link_start=1 + link_end=nhpb +! call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo + endif + write (iout,*) "molread: REFSTR",refstr if (refstr) then if (.not.pdbref) then @@ -1801,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 ! @@ -1821,6 +1924,7 @@ write(iout,*)"po setup var" subroutine pdboutC(etot,rmsd,tytul) use energy_data, only: ihpb,jhpb,itype,molnum + use energy, only:boxshift,to_box ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' @@ -1834,22 +1938,30 @@ 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) + 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) + real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3) write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),& ' ENERGY ',etot,' RMS ',rmsd iatom=0 ichain=1 ires=0 - boxxxshift(1)=int(c(1,nnt)/boxxsize) +! boxxxshift(1)=int(c(1,nnt)/boxxsize) ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1 - boxxxshift(2)=int(c(2,nnt)/boxzsize) +! boxxxshift(2)=int(c(2,nnt)/boxzsize) ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1 - boxxxshift(3)=int(c(3,nnt)/boxzsize) +! boxxxshift(3)=int(c(3,nnt)/boxzsize) ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1 - +! do i=1,3 +! if (boxxxshift(i).gt.2) boxxxshift(i)=2 +! if (boxxxshift(i).lt.-2) boxxxshift(i)=-2 +! +! enddo + do j=1,3 + cbeg(j)=c(j,nnt) + enddo + call to_box(cbeg(1),cbeg(2),cbeg(3)) boxxxx(1)=boxxsize boxxxx(2)=boxysize boxxxx(3)=boxzsize @@ -1857,25 +1969,47 @@ write(iout,*)"po setup var" do i=nnt,nct mnum=molnum(i) iti=itype(i,mnum) + if (i.ne.nct) then + iti1=itype(i+1,mnum) + mnum1=molnum(i+1) + if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle + endif if (iti.eq.ntyp1_molec(mnum)) then ichain=ichain+1 ires=0 write (ipdb,'(a)') 'TER' else ires=ires+1 +! 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 - if (mnum.eq.5) then - do j=1,3 - if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then - c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j) - else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then - c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j) - else - c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j) - endif - enddo - endif + call to_box(c(1,i),c(2,i),c(3,i)) + + DO k=1,3 + Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k)) + c(k,i)=cbeg(k)+Rdist(k) + ENDDO + if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then + call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres)) + DO k=1,3 + Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k)) + c(k,i+nres)=cbeg(k)+Rdist(k) + ENDDO + endif write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),& ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i) if ((iti.ne.10).and.(mnum.ne.5)) then