dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
- sigmasctube,krad2
- integer :: ichir1,ichir2,ijunk
+ sigmasctube,krad2,ract
+ integer :: ichir1,ichir2,ijunk,irdiff
character*3 string
- character*80 temp1
+ 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)
dsc_inv(i)=1.0D0/dsc(i)
endif
enddo
+! ip(1)=0.0001d0
+! isc(:,1)=0.0001d0
#endif
read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
do i=1,ntyp_molec(2)
- 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 *,msc(i,5),restok(i,5)
enddo
- ip(5)=0.2
+! ip(5)=0.2
! isc(5)=0.2
read (iion,*) ncatprotparm
allocate(catprm(ncatprotparm,4))
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
!
enddo
endif
enddo
+#ifdef SC_END
+ allocate(nterm_scend(2,ntyp))
+ allocate(arotam_end(0:6,2,ntyp))
+ nterm_scend=0
+ arotam_end=0.0d0
+ read (irotam_end,*) ijunk
+!c write (iout,*) "ijunk",ijunk
+ do i=1,ntyp
+ if (i.eq.10) cycle
+ do j=1,2
+ read (irotam_end,'(a)')
+ read (irotam_end,*) nterm_scend(j,i)
+!c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
+ do k=0,nterm_scend(j,i)
+ read (irotam_end,*) ijunk,arotam_end(k,j,i)
+!c write (iout,*) "k",k," arotam",arotam_end(k,j,i)
+ enddo
+ enddo
+ enddo
+!c lprint=.true.
+ if (lprint) then
+ write (iout,'(a)') &
+ "Parameters of the local potentials of sidechain ends"
+ do i=1,ntyp
+ write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
+ restyp(i,1),restyp(i,1)
+ do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
+ write (iout,'(i5,2f20.10)') &
+ j,arotam_end(j,1,i),arotam_end(j,2,i)
+ enddo
+ enddo
+ endif
+!c lprint=.false.
+#endif
+
!---------reading nucleic acid parameters for rotamers-------------------
allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
do i=1,ntyp_molec(2)
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
si=-1.0d0
do k=1,nterm_sccor(i,j)
+ print *,"test",i,j,k,l
read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
v2sccor(k,l,i,j)
v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
sigiso2(i,j)=sigiso2(j,i)
! print *,"ATU",sigma(j,i),sigma(i,j),i,j
nstate(i,j) = nstate(j,i)
- dtail(1,i,j) = dtail(1,j,i)
- dtail(2,i,j) = dtail(2,j,i)
+ dtail(1,i,j) = dtail(2,j,i)
+ dtail(2,i,j) = dtail(1,j,i)
DO k = 1, 4
alphasur(k,i,j) = alphasur(k,j,i)
wstate(k,i,j) = wstate(k,j,i)
!HERE THE MASS of MARTINI
write(*,*) "before MARTINI PARAM"
do i=1,ntyp_molec(4)
- msc(i,4)=0.0d0
- mp(4)=72.0d0
+ msc(i,4)=72.0d0
+ mp(4)=0.0d0
isc(i,4)=0.d0
enddo
ip(4)=0.0
+ msc(ntyp_molec(4)+1,4)=0.1d0
!relative dielectric constant = 15 for implicit screening
k_coulomb_lip=332.0d0/15.0d0
!kbond = 1250 kJ/(mol*nm*2)
! Ions by Aga
- allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
- allocate(alphapolcat2(ntyp,ntyp))
- allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
- allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
- allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
- allocate(epsintabcat(ntyp,ntyp))
- allocate(dtailcat(2,ntyp,ntyp))
- allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
- allocate(wqdipcat(2,ntyp,ntyp))
- allocate(wstatecat(4,ntyp,ntyp))
- allocate(dheadcat(2,2,ntyp,ntyp))
- allocate(nstatecat(ntyp,ntyp))
- allocate(debaykapcat(ntyp,ntyp))
-
- if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
- if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+ allocate(alphapolcat(ntyp,-1:ntyp_molec(5)),epsheadcat(ntyp,-1:ntyp_molec(5)),sig0headcat(ntyp,-1:ntyp_molec(5)))
+ allocate(alphapolcat2(ntyp,-1:ntyp_molec(5)))
+ allocate(sigiso1cat(ntyp,-1:ntyp_molec(5)),rborn1cat(ntyp,-1:ntyp_molec(5)),rborn2cat(ntyp,-1:ntyp_molec(5)),sigmap1cat(ntyp,-1:ntyp_molec(5)))
+ allocate(sigmap2cat(ntyp,-1:ntyp_molec(5)),sigiso2cat(ntyp,-1:ntyp_molec(5)))
+ allocate(chis1cat(ntyp,-1:ntyp_molec(5)),chis2cat(ntyp,-1:ntyp_molec(5)),wquadcat(ntyp,-1:ntyp_molec(5)),chipp1cat(ntyp,-1:ntyp_molec(5)),chipp2cat(ntyp,-1:ntyp_molec(5)))
+ allocate(epsintabcat(ntyp,-1:ntyp_molec(5)))
+ allocate(dtailcat(2,ntyp,-1:ntyp_molec(5)))
+ allocate(alphasurcat(4,ntyp,-1:ntyp_molec(5)),alphisocat(4,ntyp,-1:ntyp_molec(5)))
+ allocate(wqdipcat(2,ntyp,-1:ntyp_molec(5)))
+ allocate(wstatecat(4,ntyp,-1:ntyp_molec(5)))
+ allocate(dheadcat(2,2,ntyp,-1:ntyp_molec(5)))
+ allocate(nstatecat(ntyp,-1:ntyp_molec(5)))
+ allocate(debaykapcat(ntyp,-1:ntyp_molec(5)))
+
+ if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,-1:ntyp1))
+ if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,-1:ntyp1))
! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
- if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
- if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,-1: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
read(iion,*) ijunk
endif
print *,ntyp_molec(5)
- do i=1,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
+ ! ip(5)=0.2
! mp(5)=0.2
pstok(5)=3.0
!DIR$ NOUNROLL
- do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
+ do j=-1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
+ if (j.eq.0) cycle
do i=1,ntyp
! do j=1,ntyp_molec(5)
! write (*,*) "Im in ALAB", i, " ", j
read(iion,*) &
epscat(i,j),sigmacat(i,j), &
! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
- chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+ chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), & !6
- (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+ (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&!12
! chiscat(i,j),chiscat(j,i), &
chis1cat(i,j),chis2cat(i,j), &
- nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
- dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+ nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !19 !5 w tej lini - 1 integer pierwszy
+ dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&!23
dtailcat(1,i,j),dtailcat(2,i,j), &
- epsheadcat(i,j),sig0headcat(i,j), &
+ epsheadcat(i,j),sig0headcat(i,j), &!27
!wdipcat = w1 , w2
! rborncat(i,j),rborncat(j,i),&
rborn1cat(i,j),rborn2cat(i,j),&
- (wqdipcat(k,i,j),k=1,2), &
- alphapolcat(i,j),alphapolcat2(j,i), &
+ (wqdipcat(k,i,j),k=1,2), &!31
+ alphapolcat(i,j),alphapolcat2(j,i), &!33
(alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
END DO
END DO
- allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+ allocate(aa_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)),&
+ bb_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)))
do i=1,ntyp
- do j=1,ntyp_molec(5)
+ do j=-1,ntyp_molec(5)
+ if (j.eq.0) cycle
epsij=epscat(i,j)
rrij=sigmacat(i,j)
rrij=rrij**expon
character(len=3) :: seq,res,res2
character(len=5) :: atom
character(len=80) :: card
- real(kind=8),dimension(3,20) :: sccor
+ real(kind=8),dimension(3,40) :: sccor
integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
integer :: isugar,molecprev,firstion
character*1 :: sugar
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 '.or.(atom(1:2).eq.'ZN'))) &
- 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
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
enddo
endif
-! print *,seqalingbegin,nres
+ print *,seqalingbegin,nres
if(.not.allocated(vbld)) then
allocate(vbld(2*nres))
do i=1,2*nres
allocate(dc_norm(3,0:2*nres+2))
dc_norm(:,:)=0.d0
endif
-
+ write(iout,*) "before int_from_cart"
call int_from_cart(.true.,.false.)
call sc_loc_geom(.false.)
+ write(iout,*) "after int_from_cart"
+
+
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
enddo
+ write(iout,*) "after thetaref"
! do i=1,2*nres
! vbld_inv(i)=0.d0
! vbld(i)=0.d0
! enddo
! enddo
!
+! do i=1,2*nres
+! do j=1,3
+! chomo(j,i,k)=c(j,i)
+! enddo
+! enddo
+! write(iout,*) "after chomo"
+
if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
cou=1
write (iout,*) "symetr", symetr
do i=1,nres
- lll=lll+1
+ lll=lll+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
- if (i.gt.1) then
- if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
- chain_length=lll-1
- kkk=kkk+1
+! if (i.gt.1) then
+! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
+! chain_length=lll-1
+! kkk=kkk+1
! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
- lll=1
- endif
- endif
+! lll=1
+! endif
+! endif
do j=1,3
cref(j,i,cou)=c(j,i)
cref(j,i+nres,cou)=c(j,i+nres)
endif
enddo
enddo
- write (iout,*) chain_length
- if (chain_length.eq.0) chain_length=nres
- do j=1,3
- chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
- chain_rep(j,chain_length+nres,symetr) &
- =chain_rep(j,chain_length+nres,1)
- enddo
! diagnostic
! write (iout,*) "spraw lancuchy",chain_length,symetr
! do i=1,4
dc(j,0)=c(j,1)
enddo
- if (symetr.gt.1) then
- call permut(symetr)
- nperm=1
- do i=1,symetr
- nperm=nperm*i
- enddo
- do i=1,nperm
- write(iout,*) (tabperm(i,kkk),kkk=1,4)
- enddo
- do i=1,nperm
- cou=0
- do kkk=1,symetr
- icha=tabperm(i,kkk)
- write (iout,*) i,icha
- do lll=1,chain_length
- cou=cou+1
- if (cou.le.nres) then
- do j=1,3
- kupa=mod(lll,chain_length)
- iprzes=(kkk-1)*chain_length+lll
- if (kupa.eq.0) kupa=chain_length
- write (iout,*) "kupa", kupa
- cref(j,iprzes,i)=chain_rep(j,kupa,icha)
- cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
- enddo
- endif
- enddo
- enddo
- enddo
- endif
+! if (symetr.gt.1) then
+! call permut(symetr)
+! nperm=1
+! do i=1,symetr
+! nperm=nperm*i
+! enddo
+! do i=1,nperm
+! write(iout,*) (tabperm(i,kkk),kkk=1,4)
+! enddo
+! do i=1,nperm
+! cou=0
+! do kkk=1,symetr
+! icha=tabperm(i,kkk)
+! write (iout,*) i,icha
+! do lll=1,chain_length
+! cou=cou+1
+! if (cou.le.nres) then
+! do j=1,3
+! kupa=mod(lll,chain_length)
+! iprzes=(kkk-1)*chain_length+lll
+! if (kupa.eq.0) kupa=chain_length
+! write (iout,*) "kupa", kupa
+! cref(j,iprzes,i)=chain_rep(j,kupa,icha)
+! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
+! enddo
+! endif
+! enddo
+! enddo
+! enddo
+! endif
!-koniec robienia kopii
! diag
do kkk=1,nperm
real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
integer i
-
+ usampl=.false. ! the default value of usample should be 0
+! write(iout,*) "KURWA2",usampl
nglob_csa=0
eglob_csa=1d99
nmin_csa=0
call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
call reada(controlcard,'DRMS',drms,0.1D0)
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+ read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
+ out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
+ out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
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 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)
rest = index(controlcard,"REST").gt.0
tbf = index(controlcard,"TBF").gt.0
usampl = index(controlcard,"USAMPL").gt.0
+! write(iout,*) "KURWA",usampl
mdpdb = index(controlcard,"MDPDB").gt.0
call reada(controlcard,"T_BATH",t_bath,300.0d0)
call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
! print *,"Processor",myrank," opened file ITHEP"
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
! print *,"Processor",myrank," opened file IROTAM"
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
open (ithep,file=thetname,status='old')
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
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)
open (ithep,file=thetname,status='old',action='read')
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+ call getenv_loc('ROTPAR_END',rotname_end)
+ open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
call getenv_loc('TORDPAR',tordname)
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)
end subroutine write_stat_thread
!-----------------------------------------------------------------------------
#endif
+ subroutine readpdb_template(k)
+! Read the PDB file for read_constr_homology with read2sigma
+! and convert the peptide geometry into virtual-chain geometry.
+! implicit none
+! include 'DIMENSIONS'
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.NAMES'
+! include 'COMMON.CONTROL'
+! include 'COMMON.FRAG'
+! include 'COMMON.SETUP'
+ use compare_data, only:nhfrag,nbfrag
+ integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
+ ishift_pdb,ires_ca
+ logical lprn /.false./,fail
+ real(kind=8), dimension (3):: e1,e2,e3
+ real(kind=8) :: dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
+ character*80 card
+ real(kind=8), dimension (3,40) :: sccor
+! integer rescode
+ integer, dimension (:), allocatable :: iterter
+ if(.not.allocated(iterter))allocate(iterter(nres))
+ do i=1,nres
+ iterter(i)=0
+ enddo
+ ibeg=1
+ ishift1=0
+ ishift=0
+ write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
+ do
+ read (ipdbin,'(a80)',end=10) card
+ if (card(:3).eq.'END') then
+ goto 10
+ else if (card(:3).eq.'TER') then
+! End current chain
+ ires_old=ires+2
+ itype(ires_old-1,1)=ntyp1
+ iterter(ires_old-1)=1
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ if (ires_old.lt.nres) then
+#endif
+ itype(ires_old,1)=ntyp1
+ iterter(ires_old)=1
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ endif
+#endif
+ ibeg=2
+! write (iout,*) "Chain ended",ires,ishift,ires_old
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+! Fish out the ATOM cards.
+ if (index(card(1:4),'ATOM').gt.0) then
+ read (card(12:16),*) atom
+! write (iout,*) "! ",atom," !",ires
+! if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+! write (iout,*) "ires",ires,ires-ishift+ishift1,
+! & " ires_old",ires_old
+! write (iout,*) "ishift",ishift," ishift1",ishift1
+! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
+! Calculate the CM of the preceding residue.
+ if (ibeg.eq.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires_old)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
+ endif
+! Start new residue.
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+! write (iout,*) "BEG ires",ires
+ ishift=ires-1
+ if (res.ne.'GLY' .and. res.ne. 'ACE') then
+ ishift=ishift-1
+ itype(1,1)=ntyp1
+ endif
+ ires=ires-ishift+ishift1
+ ires_old=ires
+! write (iout,*) "ishift",ishift," ires",ires,
+! & " ires_old",ires_old
+! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ibeg=0
+ else if (ibeg.eq.2) then
+! Start a new chain
+ ishift=-ires_old+ires-1
+ ires=ires_old+1
+! write (iout,*) "New chain started",ires,ishift
+ ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
+ endif
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires,1)=10
+ else
+ itype(ires,1)=rescode(ires,res,0,1)
+ endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+! write (iout,*) "ires_old",ires_old," ires",ires
+! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+! ishift1=ishift1+1
+! endif
+! write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
+ res.eq.'NHE'.and.atom(:2).eq.'HN') then
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)') &
+ ires,itype(ires,1),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
+ do j=1,3
+ sccor(j,iii)=c(j,ires)
+ enddo
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+! write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
+ atom.ne.'N' .and. atom.ne.'C' .and.&
+ atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
+ atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+! write (iout,*) "sidechain ",atom
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+ endif
+ enddo
+ 10 if(me.eq.king.or..not.out1file) &
+ write (iout,'(a,i5)') ' Nres: ',ires
+! Calculate dummy residue coordinates inside the "chain" of a multichain
+! system
+ nres=ires
+ write(2,*) "tutaj",ires,nres
+ do i=2,nres-1
+! write (iout,*) i,itype(i),itype(i+1)
+ if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
+ if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif !fail
+ do j=1,3
+ c(j,i)=c(j,i-1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i-2)-c(j,i-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i-1)+dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ else !itype(i+1).eq.ntyp1
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i+3)-c(j,i+2))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i+1)-dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ endif !itype(i+1).eq.ntyp1
+ endif !itype.eq.ntyp1
+ enddo
+! Calculate the CM of the last side chain.
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ nsup=nres
+ nstart_sup=1
+ if (itype(nres,1).ne.10) then
+ nres=nres+1
+ itype(nres,1)=ntyp1
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ endif
+ do i=2,nres-1
+ do j=1,3
+ c(j,i+nres)=dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ if (itype(1,1).eq.ntyp1) then
+ nsup=nsup-1
+ nstart_sup=2
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,4)-c(j,3))/2.0
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ endif
+! Copy the coordinates to reference coordinates
+! do i=1,2*nres
+! do j=1,3
+! cref(j,i)=c(j,i)
+! enddo
+! enddo
+! Calculate internal coordinates.
+ if (out_template_coord) then
+ write (iout,'(/a)') &
+ "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
+ "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
+ restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
+ (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+! Calculate internal coordinates.
+ call int_from_cart(.true.,out_template_coord)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+! & vbld_inv(i+nres)
+ enddo
+ do i=1,nres
+ do j=1,3
+ cref(j,i,1)=c(j,i)
+ cref(j,i+nres,1)=c(j,i+nres)
+ enddo
+ enddo
+ do i=1,2*nres
+ do j=1,3
+ chomo(j,i,k)=c(j,i)
+ enddo
+ enddo
+
+ return
+ end subroutine readpdb_template
+!-----------------------------------------------------------------------------
+!#endif
!-----------------------------------------------------------------------------
end module io_config
projects / unres4.git / blobdiff