X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fio.F90;h=b2ebc22463690f46358d9915555f12ebe0683752;hb=6838b92b256fbeb9953f1e77eed0a0ec68b12192;hp=cc8b68171acdced847af34b8351a5ed3dfedab82;hpb=df2469d9ac903d93889867f4e50e9bf6c428c1c6;p=unres4.git diff --git a/source/unres/io.F90 b/source/unres/io.F90 index cc8b681..b2ebc22 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -736,10 +736,12 @@ real(kind=8),dimension(3,maxres2+2) :: c_alloc real(kind=8),dimension(3,0:maxres2) :: dc_alloc + real(kind=8),dimension(:,:), allocatable :: secprob integer,dimension(maxres) :: itype_alloc integer :: iti,nsi,maxsi,itrial,itmp - real(kind=8) :: wlong,scalscp,co,ssscale + real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,& + sigmabet,sumv allocate(weights(n_ene)) !----------------------------- allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres @@ -782,9 +784,10 @@ 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) +! print *,"KUR..",wtor_nucl 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,'WCORR3_NUC',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) @@ -830,6 +833,26 @@ weights(17)=wbond weights(18)=scal14 weights(21)=wsccor + weights(26)=wvdwpp_nucl + weights(27)=welpp + weights(28)=wvdwpsb + weights(29)=welpsb + weights(30)=wvdwsb + weights(31)=welsb + weights(32)=wbond_nucl + weights(33)=wang_nucl + weights(34)=wsbloc + weights(35)=wtor_nucl + weights(36)=wtor_d_nucl + weights(37)=wcorr_nucl + weights(38)=wcorr3_nucl + weights(41)=wcatcat + weights(42)=wcatprot + weights(46)=wscbase + weights(47)=wpepbase + weights(48)=wscpho + weights(49)=wpeppho + if(me.eq.king.or..not.out1file) & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,& @@ -953,6 +976,8 @@ 34 continue ! print *,'Begin reading pdb data' call readpdb +! call int_from_cart1(.true.) + ! print *,'Finished reading pdb data' if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup,& @@ -1031,7 +1056,7 @@ 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) + itype(itmp,molec)=rescode(i,sequence(i,molec)(1:2),iscode,molec) enddo endif @@ -1052,7 +1077,7 @@ nres=0 do i=1,5 nres=nres+nres_molec(i) - print *,nres_molec(i) + print *,"nres_molec",nres,nres_molec(i) enddo ! Assign initial virtual bond lengths @@ -1067,14 +1092,18 @@ enddo ! print *,nres_molec(i) endif - if(.not.allocated(vbld)) allocate(vbld(2*nres)) + print *,nres,"nres" + if(.not.allocated(vbld)) then + print *, "I DO ENTER" + allocate(vbld(2*nres)) + endif if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres)) do i=2,nres vbld(i)=vbl vbld_inv(i)=vblinv enddo do i=2,nres-1 - print *, "molnum",molnum(i),itype(i,molnum(i)) +! print *, "molnum",molnum(i),itype(i,molnum(i)),nres,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), @@ -1143,7 +1172,7 @@ endif call read_bridge !-------------------------------- - print *,"tu dochodze" +! print *,"tu dochodze" ! znamy nres oraz nss można zaalokowac potrzebne tablice call alloc_geo_arrays call alloc_ener_arrays @@ -1155,6 +1184,7 @@ enddo read (inp,*) ndih_constr if (ndih_constr.gt.0) then + raw_psipred=.false. allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr) allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr) allocate(ftors(ndih_constr)) !(maxdih_constr) @@ -1181,6 +1211,45 @@ phibound(1,ii) = phi0(i)-drange(i) phibound(2,ii) = phi0(i)+drange(i) enddo + else if (ndih_constr.lt.0) then + raw_psipred=.true. + allocate(secprob(3,nres)) + allocate(vpsipred(3,nres)) + allocate(sdihed(2,nres)) + call card_concat(weightcard,.true.) + call reada(weightcard,"PHIHEL",phihel,50.0D0) + call reada(weightcard,"PHIBET",phibet,180.0D0) + call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0) + call reada(weightcard,"SIGMABET",sigmabet,40.0d0) + call reada(weightcard,"WDIHC",wdihc,0.591D0) + write (iout,*) "Weight of dihedral angle restraints",wdihc + read(inp,'(9x,3f7.3)') & + (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) + write (iout,*) "The secprob array" + do i=nnt,nct + write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3) + enddo + ndih_constr=0 + do i=nnt+3,nct + if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 & + .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then + ndih_constr=ndih_constr+1 + idih_constr(ndih_constr)=i + sumv=0.0d0 + do j=1,3 + vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2) + sumv=sumv+vpsipred(j,ndih_constr) + enddo + do j=1,3 + vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv + enddo + phibound(1,ndih_constr)=phihel*deg2rad + phibound(2,ndih_constr)=phibet*deg2rad + sdihed(1,ndih_constr)=sigmahel*deg2rad + sdihed(2,ndih_constr)=sigmabet*deg2rad + endif + enddo + endif if (with_theta_constr) then !C with_theta_constr is keyword allowing for occurance of theta constrains @@ -1371,7 +1440,10 @@ enddo return else - call read_angles(inp,*36) + write(iout,*) "read angles from input" + call read_angles(inp,*36) + call chainbuild + endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' @@ -1508,6 +1580,10 @@ enddo endif if (i2ndstr.gt.0) call secstrp2dihc + if (indpdb.gt.0) then + write(iout,*) "WCHODZE TU!!" + call int_from_cart1(.true.) + endif ! call geom_to_var(nvar,x) ! call etotal(energia(0)) ! call enerprint(energia(0))