include 'COMMON.CONTROL'
include 'COMMON.SBRIDGE'
include 'COMMON.IOUNITS'
+ include 'COMMON.SPLITELE'
logical file_exist
C Read force-field parameters except weights
call parmread
include 'COMMON.FFIELD'
include 'COMMON.INTERACT'
include 'COMMON.SETUP'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SHIELD'
COMMON /MACHSW/ KDIAG,ICORFL,IXDR
character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
character*80 ucase
C Set up the time limit (caution! The time must be input in minutes!)
read_cart=index(controlcard,'READ_CART').gt.0
call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+C this variable with_theta_constr is the variable which allow to read and execute the
+C constrains on theta angles WITH_THETA_CONSTR is the keyword
+ with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+ write (iout,*) "with_theta_constr ",with_theta_constr
call readi(controlcard,'SYM',symetr,1)
- call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
+ call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
indpdb=index(controlcard,'PDBSTART')
extconf=(index(controlcard,'EXTCONF').gt.0)
+ AFMlog=(index(controlcard,'AFM'))
+ selfguide=(index(controlcard,'SELFGUIDE'))
+ print *,'AFMlog',AFMlog,selfguide,"KUPA"
+ call readi(controlcard,'TUBEMOD',tubelog,0)
+ write (iout,*) TUBElog,"TUBEMODE"
+ call readi(controlcard,'GENCONSTR',genconstr,0)
+C write (iout,*) TUBElog,"TUBEMODE"
call readi(controlcard,'IPRINT',iprint,0)
+C SHIELD keyword sets if the shielding effect of side-chains is used
+C 0 denots no shielding is used all peptide are equally despite the
+C solvent accesible area
+C 1 the newly introduced function
+C 2 reseved for further possible developement
+ call readi(controlcard,'SHIELD',shield_mode,0)
+C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "shield_mode",shield_mode
+C endif
+ 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 readi(controlcard,'MAXGEN',maxgen,10000)
call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
call readi(controlcard,"KDIAG",kdiag,0)
else
icheckgrad=3
endif
+ call reada(controlcard,'DELTA',aincr,1.0d-7)
elseif (index(controlcard,'THREAD').gt.0) then
modecalc=2
call readi(controlcard,'THREAD',nthread,0)
i2ndstr=index(controlcard,'USE_SEC_PRED')
gradout=index(controlcard,'GRADOUT').gt.0
gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
+C DISTCHAINMAX become obsolete for periodic boundry condition
call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
+C Reading the dimensions of box in x,y,z coordinates
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+c Cutoff range for interactions
+ call reada(controlcard,"R_CUT",r_cut,15.0d0)
+ call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ 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
+C endif
+ 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
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+ write (iout,*) "SHIELD MODE",shield_mode
+ if (TUBElog.gt.0) then
+ call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+ call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+ call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+ call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
+ call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
+ call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
+ buftubebot=bordtubebot+tubebufthick
+ buftubetop=bordtubetop-tubebufthick
+ endif
+ if (shield_mode.gt.0) then
+ pi=3.141592d0
+C VSolvSphere the volume of solving sphere
+C print *,pi,"pi"
+C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
+C there will be no distinction between proline peptide group and normal peptide
+C group in case of shielding parameters
+ VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+ VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+ write (iout,*) VSolvSphere,VSolvSphere_div
+C long axis of side chain
+ do i=1,ntyp
+ long_r_sidechain(i)=vbldsc0(1,i)
+ short_r_sidechain(i)=sigma0(i)
+ enddo
+ buff_shield=1.0d0
+ endif
if (me.eq.king .or. .not.out1file )
& write (iout,*) "DISTCHAINMAX",distchainmax
ntime_split0=ntime_split
call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
ntime_split0=ntime_split
- call reada(controlcard,"R_CUT",r_cut,2.0d0)
- call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+c call reada(controlcard,"R_CUT",r_cut,2.0d0)
+c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
rest = index(controlcard,"REST").gt.0
tbf = index(controlcard,"TBF").gt.0
usampl = index(controlcard,"USAMPL").gt.0
include 'COMMON.BOUNDS'
include 'COMMON.MD'
include 'COMMON.SETUP'
+ include 'COMMON.SHIELD'
character*4 sequence(maxres)
integer rescode
double precision x(maxvar)
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,'WSHIELD',wshield,1.0d0)
+ call reada(weightcard,'WLT',wliptran,0.0D0)
+ call reada(weightcard,'WTUBE',wtube,1.0D0)
if (index(weightcard,'SOFT').gt.0) ipot=6
C 12/1/95 Added weight for the multi-body term WCORR
call reada(weightcard,'WCORRH',wcorr,1.0D0)
call reada(weightcard,"BTRISS",btriss,0.021D0)
call reada(weightcard,"CTRISS",ctriss,1.001D0)
call reada(weightcard,"DTRISS",dtriss,1.001D0)
+ call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
write (iout,*) "ATRISS=", atriss
write (iout,*) "BTRISS=", btriss
write (iout,*) "CTRISS=", ctriss
write (iout,*) "DTRISS=", dtriss
+ if (shield_mode.gt.0) then
+ lipscale=0.0d0
+ write (iout,*) "liscale not used in shield mode"
+ endif
+ write (iout,*) "lipid scaling factor", lipscale
dyn_ss=(index(weightcard,'DYN_SS').gt.0)
do i=1,maxres
dyn_ss_mask(i)=.false.
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
+ akcm=akcm/wsc
+ akth=akth/wsc
+ akct=akct/wsc
+ v1ss=v1ss/wsc
+ v2ss=v2ss/wsc
+ v3ss=v3ss/wsc
else
if (wstrain.ne.0.0) then
ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
ss_depth=0.0
endif
endif
-
+ write (iout,*) "wshield,", wshield
if(me.eq.king.or..not.out1file) then
write (iout,*) "Parameters of the SS-bond potential:"
write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
phibound(2,ii) = phi0(i)+drange(i)
enddo
endif
+C first setting the theta boundaries to 0 to pi
+C this mean that there is no energy penalty for any angle occuring this can be applied
+C for generate random conformation but is not implemented in this way
+C do i=1,nres
+C thetabound(1,i)=0
+C thetabound(2,i)=pi
+C enddo
+C begin reading theta constrains this is quartic constrains allowing to
+C have smooth second derivative
+ 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
+ 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
+C
+C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
+C write (iout,*) "with_dihed_constr ",with_dihed_constr
nnt=1
#ifdef MPI
if (me.eq.king) then
call MPI_Finalize(MPI_COMM_WORLD,IERROR)
stop 'Error reading reference structure'
#endif
- 39 call chainbuild
+ 39 call chainbuild_extconf
call setup_var
czscore call geom_to_var(nvar,coord_exp_zs(1,1))
nstart_sup=nnt
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
return
else
call read_angles(inp,*36)
+ call chainbuild_extconf
endif
goto 37
36 write (iout,'(a)') 'Error reading angle file.'
omeg(i)=-120d0*deg2rad
if (itype(i).le.0) omeg(i)=-omeg(i)
enddo
+ call chainbuild_extconf
else
if(me.eq.king.or..not.out1file)
& write (iout,'(a)') 'Random-generated initial geometry.'
40 continue
endif
#else
- do itrial=1,100
- itmp=1
- call gen_rand_conf(itmp,*30)
- goto 40
- 30 write (iout,*) 'Failed to generate random conformation',
- & ', itrial=',itrial
- write (*,*) 'Failed to generate random conformation',
- & ', itrial=',itrial
- enddo
- write (iout,'(a,i3,a)') 'Processor:',me,
- & ' error in generating random conformation.'
- write (*,'(a,i3,a)') 'Processor:',me,
+ write (*,'(a)')
& ' error in generating random conformation.'
stop
40 continue
cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
cd & ' nsup',nsup
+ if (constr_dist.eq.11) then
+ do i=nstart_sup,nstart_sup+nsup-1
+ do j=i+2,nstart_sup+nsup-1
+ distance=dist(i,j)
+ if (distance.le.15.0) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=i+nstart_seq-nstart_sup
+ jhpb(nhpb)=j+nstart_seq-nstart_sup
+ forcon(nhpb)=sqrt(0.04*distance)
+ fordepth(nhpb)=sqrt(40.0/distance)
+ dhpb(nhpb)=distance-0.1d0
+ dhpb1(nhpb)=distance+0.1d0
+
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+ write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+ endif
+ enddo
+ enddo
+ else
do i=nstart_sup,nstart_sup+nsup-1
cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
cd & ' seq_pdb', restyp(itype_pdb(i))
forcon(nhpb)=weidis
dhpb(nhpb)=dist(i,j)
enddo
- enddo
+ enddo
+ endif
cd write (iout,'(a)') 'Distance constraints:'
cd do i=nss+1,nhpb
cd ii=ihpb(i)
open (itorp,file=torname,status='old',readonly,shared)
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old',readonly,shared)
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old',readonly,shared)
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old',readonly,shared)
call getenv_loc('FOURIER',fouriername)
open (ifourier,file=fouriername,status='old',readonly,shared)
call getenv_loc('ELEPAR',elename)
open (ielep,file=elename,status='old',readonly,shared)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly,shared)
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',readonly,shared)
#elif (defined CRAY) || (defined AIX)
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
& action='read')
c print *,"Processor",myrank," opened file ITORP"
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old',action='read')
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old',action='read')
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old',action='read')
c print *,"Processor",myrank," opened file ITORDP"
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old',action='read')
c print *,"Processor",myrank," opened file IELEP"
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',action='read')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
c print *,"Processor",myrank," opened file ISIDEP"
c print *,"Processor",myrank," opened parameter files"
#elif (defined G77)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old')
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old')
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old')
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old')
call getenv_loc('FOURIER',fouriername)
open (ielep,file=elename,status='old')
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old')
#else
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
& readonly)
open (itorp,file=torname,status='old',readonly)
call getenv_loc('TORDPAR',tordname)
open (itordp,file=tordname,status='old',readonly)
+ call getenv_loc('TORKCC',torkccname)
+ open (itorkcc,file=torkccname,status='old',readonly)
+ call getenv_loc('THETKCC',thetkccname)
+ open (ithetkcc,file=thetkccname,status='old',readonly)
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old',readonly)
#ifndef CRYST_THETA
open (ielep,file=elename,status='old',readonly)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly)
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
#ifndef CRYST_SC
call getenv_loc('ROTPARPDB',rotname_pdb)
open (irotam_pdb,file=rotname_pdb,status='old',action='read')
#endif
#endif
+ call getenv_loc('TUBEPAR',tubename)
+ open (itube,file=tubename,status='old',readonly)
#ifndef OLDSCP
C
C 8/9/01 In the newest version SCp interaction constants are read from a file
mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
if (lentmp.gt.0)
- & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
& //'.stat')
rest2name=prefix(:ilen(prefix))//'.rst'
if(usampl) then
include 'COMMON.MD'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
+ totTafm=totT
do i=1,2*nres
read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
enddo
return
end
+C---------------------------------------------------------------------------
+ subroutine read_afminp
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ character*320 afmcard
+ print *, "wchodze"
+ call card_concat(afmcard)
+ call readi(afmcard,"BEG",afmbeg,0)
+ call readi(afmcard,"END",afmend,0)
+ call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
+ call reada(afmcard,"VEL",velAFMconst,0.0d0)
+ print *,'FORCE=' ,forceAFMconst
+CCCC NOW PROPERTIES FOR AFM
+ distafminit=0.0d0
+ do i=1,3
+ distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
+ enddo
+ distafminit=dsqrt(distafminit)
+ print *,'initdist',distafminit
+ return
+ end
+
c-------------------------------------------------------------------------------
subroutine read_dist_constr
implicit real*8 (a-h,o-z)
write (iout,*) "Calling read_dist_constr"
c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
c call flush(iout)
+ if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
+ call gen_dist_constr
+ go to 1712
+ endif
call card_concat(controlcard)
call readi(controlcard,"NFRAG",nfrag_,0)
call readi(controlcard,"NPAIR",npair_,0)
#endif
enddo
+ 1712 continue
call flush(iout)
return
end
projects / unres.git / blobdiff