X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio.f90;h=ab5dc2aef24fbfaab3c9fc1151ac8a18f99c594c;hb=7077a079d2932fdbb6cb5c28bae95bad9ff9105c;hp=7964cdd3c53a593df9e699a9b1d698489cf9073c;hpb=5d299c1a16ab51f8206b8ee3b17c7bcabe9321b7;p=unres4.git diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 7964cdd..ab5dc2a 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -220,12 +220,12 @@ DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL) do i=nnt,nct - if (itype(i).ne.10) then -!d print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1) + if (itype(i,1).ne.10) then +!d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1) fail=.true. ii=0 do while (fail .and. ii .le. maxsi) - call gen_side(itype(i),theta(i+1),alph(i),omeg(i),fail) + call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail) ii = ii+1 enddo endif @@ -487,7 +487,7 @@ subroutine statout(itime) use energy_data - use control_data, only:refstr + use control_data use MD_data use MPI_data use compare, only:rms_nac_nnc @@ -513,8 +513,8 @@ character(len=4) :: format1,format2 character(len=30) :: format !el local variables - integer :: i,ii1,ii2 - real(kind=8) :: rms,frac,frac_nn,co + integer :: i,ii1,ii2,j + real(kind=8) :: rms,frac,frac_nn,co,distance #ifdef AIX if(itime.eq.0) then @@ -527,6 +527,46 @@ open(istat,file=statname,access="append") #endif #endif + if (AFMlog.gt.0) then + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')& + itime,totT,EK,potE,totE,& + rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& + potEcomp(23),me + format1="a133" + else +!C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,7f12.3,i5,$)') & + itime,totT,EK,potE,totE,& + kinetic_T,t_bath,gyrate(),& + potEcomp(23),me + format1="a114" + endif + else if (selfguide.gt.0) then + distance=0.0 + do j=1,3 + distance=distance+(c(j,afmend)-c(j,afmbeg))**2 + enddo + distance=dsqrt(distance) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, & + f9.3,i5,$)') & + itime,totT,EK,potE,totE,& + rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),& + distance,potEcomp(23),me + format1="a133" +!C print *,"CHUJOWO" + else +!C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,8f12.3,i5,$)')& + itime,totT,EK,potE,totE, & + kinetic_T,t_bath,gyrate(),& + distance,potEcomp(23),me + format1="a114" + endif + else if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') & @@ -539,6 +579,7 @@ amax,kinetic_T,t_bath,gyrate(),me format1="a114" endif + ENDIF if(usampl.and.totT.gt.eq_time) then write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,& (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),& @@ -675,11 +716,11 @@ ! include 'COMMON.BOUNDS' ! include 'COMMON.MD' ! include 'COMMON.SETUP' - character(len=4),dimension(:),allocatable :: sequence !(maxres) + character(len=4),dimension(:,:),allocatable :: sequence !(maxres,maxmolecules) ! integer :: rescode ! double precision x(maxvar) character(len=256) :: pdbfile - character(len=320) :: weightcard + character(len=800) :: weightcard character(len=80) :: weightcard_t!,ucase ! integer,dimension(:),allocatable :: itype_pdb !(maxres) ! common /pizda/ itype_pdb @@ -695,18 +736,19 @@ integer,dimension(maxres) :: itype_alloc integer :: iti,nsi,maxsi,itrial,itmp - real(kind=8) :: wlong,scalscp,co + real(kind=8) :: wlong,scalscp,co,ssscale allocate(weights(n_ene)) !----------------------------- allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres allocate(dc(3,0:2*maxres)) !(3,0:maxres2) - allocate(itype(maxres)) !(maxres) + allocate(itype(maxres,5)) !(maxres) + allocate(istype(maxres)) ! ! Zero out tables. ! c(:,:)=0.0D0 dc(:,:)=0.0D0 - itype(:)=0 + itype(:,:)=0 !----------------------------- ! ! Body @@ -727,12 +769,26 @@ call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) + call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0) + call reada(weightcard,'WELPP',welpp,0.0d0) + call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0) + call reada(weightcard,'WELPSB',welpsb,0.0D0) + call reada(weightcard,'WVDWSB',wvdwsb,0.0d0) + call reada(weightcard,'WELSB',welsb,0.0D0) + call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0) + call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0) + call reada(weightcard,'WSBLOC',wsbloc,0.0D0) + call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0) + call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0) + call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0) + call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,0.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WSHIELD',wshield,0.05D0) call reada(weightcard,'LIPSCALE',lipscale,1.0D0) call reada(weightcard,'WLT',wliptran,1.0D0) + call reada(weightcard,'WTUBE',wtube,1.0d0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) @@ -740,7 +796,15 @@ call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) + call reada(weightcard,'WSCBASE',wscbase,0.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 + call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0) + call reada(weightcard,'WCATCAT',wcatcat,0.0d0) + call reada(weightcard,'WCATPROT',wcatprot,0.0d0) + call reada(weightcard,'WPEPBASE',wpepbase,1.0d0) + call reada(weightcard,'WSCPHO',wscpho,0.0d0) + call reada(weightcard,'WPEPPHO',wpeppho,0.0d0) + ! 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 @@ -844,18 +908,23 @@ 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,"HT",Ht,0.0D0) if (dyn_ss) then - ss_depth=ebr/wsc-0.25*eps(1,1) - Ht=Ht/wsc-0.25*eps(1,1) - akcm=akcm*wstrain/wsc - akth=akth*wstrain/wsc - akct=akct*wstrain/wsc - v1ss=v1ss*wstrain/wsc - v2ss=v2ss*wstrain/wsc - v3ss=v3ss*wstrain/wsc + 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 @@ -888,13 +957,13 @@ !el if(.not.allocated(itype_pdb)) allocate(itype_pdb(nres)) do i=1,nres - itype_pdb(i)=itype(i) + itype_pdb(i)=itype(i,1) enddo close (ipdbin) nnt=nstart_sup nct=nstart_sup+nsup-1 !el if(.not.allocated(icont_ref)) - allocate(icont_ref(2,12*nres)) ! maxcont=12*maxres + allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then @@ -902,8 +971,8 @@ write(iout,*)'Adding sidechains' maxsi=1000 do i=2,nres-1 - iti=itype(i) - if (iti.ne.10 .and. itype(i).ne.ntyp1) then + iti=itype(i,1) + if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) @@ -916,40 +985,101 @@ enddo endif endif - + if (indpdb.eq.0) then + nres_molec(:)=0 + allocate(sequence(maxres,5)) +! itype(:,:)=0 + itmp=0 + if (protein) then ! Read sequence if not taken from the pdb file. - read (inp,*) nres -! print *,'nres=',nres - allocate(sequence(nres)) + molec=1 + read (inp,*) nres_molec(molec) + print *,'nres=',nres if (iscode.gt.0) then - read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) + read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec)) else - read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) + read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec)) endif +! read(inp,*) weightcard_t +! print *,"po seq" weightcard_t ! Convert sequence to numeric code - do i=1,nres - itype(i)=rescode(i,sequence(i),iscode) + + do i=1,nres_molec(molec) + itmp=itmp+1 + itype(i,1)=rescode(i,sequence(i,molec),iscode,molec) + print *,itype(i,1) + + enddo + endif +! read(inp,*) weightcard_t +! print *,"po seq", weightcard_t + + if (nucleic) then +! Read sequence if not taken from the pdb file. + molec=2 + read (inp,*) nres_molec(molec) +! print *,'nres=',nres +! allocate(sequence(maxres,5)) +! if (iscode.gt.0) then + read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec)) +! Convert sequence to numeric code + + do i=1,nres_molec(molec) + itmp=itmp+1 + istype(itmp)=sugarcode(sequence(i,molec)(1:1),i) + itype(itmp,molec)=rescode(i,sequence(i,molec)(2:2),iscode,molec) + enddo + endif + + if (ions) then +! Read sequence if not taken from the pdb file. + molec=5 + read (inp,*) nres_molec(molec) +! print *,'nres=',nres + read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec)) +! Convert sequence to numeric code + print *,nres_molec(molec) + do i=1,nres_molec(molec) + itmp=itmp+1 + print *,itmp,"itmp" + itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec) enddo + endif + nres=0 + do i=1,5 + nres=nres+nres_molec(i) + print *,nres_molec(i) + enddo + ! Assign initial virtual bond lengths -!elwrite(iout,*) "test_alloc" + if(.not.allocated(molnum)) then + allocate(molnum(nres+1)) + itmp=0 + do i=1,5 + do j=1,nres_molec(i) + itmp=itmp+1 + molnum(itmp)=i + enddo + enddo +! print *,nres_molec(i) + endif if(.not.allocated(vbld)) allocate(vbld(2*nres)) -!elwrite(iout,*) "test_alloc" if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres)) -!elwrite(iout,*) "test_alloc" do i=2,nres vbld(i)=vbl vbld_inv(i)=vblinv enddo do i=2,nres-1 - vbld(i+nres)=dsc(iabs(itype(i))) - vbld_inv(i+nres)=dsc_inv(iabs(itype(i))) -! write (iout,*) "i",i," itype",itype(i), -! & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres) + print *, "molnum",molnum(i),itype(i,molnum(i)) + vbld(i+nres)=dsc(iabs(itype(i,molnum(i)))) + vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i)))) +! write (iout,*) "i",i," itype",itype(i,1), +! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres) enddo endif ! print *,nres -! print '(20i4)',(itype(i),i=1,nres) +! print '(20i4)',(itype(i,1),i=1,nres) !---------------------------- !el reallocate tables ! do i=1,maxres2 @@ -959,8 +1089,8 @@ ! enddo ! enddo ! do i=1,maxres -!elwrite(iout,*) "itype",i,itype(i) -! itype_alloc(i)=itype(i) +!elwrite(iout,*) "itype",i,itype(i,1) +! itype_alloc(i)=itype(i,1) ! enddo ! deallocate(c) @@ -979,21 +1109,21 @@ ! enddo ! enddo ! do i=1,nres+2 -! itype(i)=itype_alloc(i) +! itype(i,1)=itype_alloc(i) ! itel(i)=0 ! enddo !-------------------------- do i=1,nres #ifdef PROCOR - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then #else - if (itype(i).eq.ntyp1) then + if (itype(i,1).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR - else if (iabs(itype(i+1)).ne.20) then + else if (iabs(itype(i+1,1)).ne.20) then #else - else if (iabs(itype(i)).ne.20) then + else if (iabs(itype(i,1)).ne.20) then #endif itel(i)=1 else @@ -1002,13 +1132,15 @@ enddo if(me.eq.king.or..not.out1file)then write (iout,*) "ITEL" + print *,nres,"nres" do i=1,nres-1 - write (iout,*) i,itype(i),itel(i) + write (iout,*) i,itype(i,1),itel(i) enddo print *,'Call Read_Bridge.' endif call read_bridge !-------------------------------- + print *,"tu dochodze" ! znamy nres oraz nss można zaalokowac potrzebne tablice call alloc_geo_arrays call alloc_ener_arrays @@ -1022,6 +1154,7 @@ if (ndih_constr.gt.0) then allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr) allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr) + read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) if(me.eq.king.or..not.out1file)then @@ -1043,6 +1176,47 @@ phibound(2,ii) = phi0(i)+drange(i) enddo endif + if (with_theta_constr) then +!C with_theta_constr is keyword allowing for occurance of theta constrains + read (inp,*) ntheta_constr +!C ntheta_constr is the number of theta constrains + if (ntheta_constr.gt.0) then +!C read (inp,*) ftors + allocate(itheta_constr(ntheta_constr)) + allocate(theta_constr0(ntheta_constr)) + allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr)) + read (inp,*) (itheta_constr(i),theta_constr0(i), & + theta_drange(i),for_thet_constr(i), & + i=1,ntheta_constr) +!C the above code reads from 1 to ntheta_constr +!C itheta_constr(i) residue i for which is theta_constr +!C theta_constr0 the global minimum value +!C theta_drange is range for which there is no energy penalty +!C for_thet_constr is the force constant for quartic energy penalty +!C E=k*x**4 + if(me.eq.king.or..not.out1file)then + write (iout,*) & + 'There are',ntheta_constr,' constraints on phi angles.' + do i=1,ntheta_constr + write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), & + theta_drange(i), & + for_thet_constr(i) + enddo + endif + do i=1,ntheta_constr + theta_constr0(i)=deg2rad*theta_constr0(i) + theta_drange(i)=deg2rad*theta_drange(i) + enddo +!C if(me.eq.king.or..not.out1file) +!C & write (iout,*) 'FTORS',ftors +!C do i=1,ntheta_constr +!C ii = itheta_constr(i) +!C thetabound(1,ii) = phi0(i)-drange(i) +!C thetabound(2,ii) = phi0(i)+drange(i) +!C enddo + endif ! ntheta_constr.gt.0 + endif! with_theta_constr + nnt=1 #ifdef MPI if (me.eq.king) then @@ -1050,22 +1224,22 @@ write (iout,'(a)') 'Boundaries in phi angle sampling:' do i=1,nres write (iout,'(a3,i5,2f10.1)') & - restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg + restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg enddo #ifdef MP endif #endif nct=nres -!d print *,'NNT=',NNT,' NCT=',NCT - if (itype(1).eq.ntyp1) nnt=2 - if (itype(nres).eq.ntyp1) nct=nct-1 + print *,'NNT=',NNT,' NCT=',NCT + if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2 + if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1 if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup nstart_seq=nnt if (nsup.le.(nct-nnt+1)) then do i=0,nct-nnt+1-nsup - if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then + if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then nstart_seq=nnt+i goto 111 endif @@ -1075,7 +1249,7 @@ stop else do i=0,nsup-(nct-nnt+1) - if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) & + if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) & then nstart_sup=nstart_sup+i nsup=nct-nnt+1 @@ -1133,6 +1307,7 @@ call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb + if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co @@ -1144,9 +1319,9 @@ icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup enddo if(me.eq.king.or..not.out1file) & - write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',& + write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',& icont_ref(1,i),' ',& - restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) + restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i) enddo endif endif @@ -1164,7 +1339,8 @@ ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "Exit READ_CART" write (iout,'(8f10.5)') & - ((c(l,k),l=1,3),k=1,nres),& + ((c(l,k),l=1,3),k=1,nres) + write (iout,'(8f10.5)') & ((c(l,k+nres),l=1,3),k=nnt,nct) call int_from_cart1(.true.) write (iout,*) "Finish INT_TO_CART" @@ -1175,7 +1351,7 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres) @@ -1205,10 +1381,9 @@ do i=2,nres-1 alph(i)=110d0*deg2rad enddo -!elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres) do i=2,nres-1 omeg(i)=-120d0*deg2rad - if (itype(i).le.0) omeg(i)=-omeg(i) + if (itype(i,1).le.0) omeg(i)=-omeg(i) enddo else if(me.eq.king.or..not.out1file) & @@ -1283,10 +1458,8 @@ call read_threadbase endif call setup_var -!elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres) if (me.eq.king .or. .not. out1file) & call intout -!elwrite (iout,*)"alph(i)*rad2deg",(alph(i)*rad2deg, i=1,nres) if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' @@ -1298,11 +1471,11 @@ do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres - it1=itype(i1) - it2=itype(i2) + it1=itype(i1,1) + it2=itype(i2,1) if (me.eq.king.or..not.out1file) & write (iout,'(2a,i3,3a,i3,a,3f10.3)') & - restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),& + restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),& ebr,forcon(i) enddo write (iout,'(a)')