X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=89ddeaca181c1bb74aa8bbc8bf291f8661211f37;hb=0cda1be1a3786bb92473b33bf4694c86e35e2407;hp=286cc57cf99b62749c39b36fd5e16aded298d624;hpb=31c90fb67d7b1da11415fe4cb5559320182d00d8;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 286cc57..89ddeac 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -8,6 +8,7 @@ include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' + include 'COMMON.SPLITELE' logical file_exist C Read force-field parameters except weights call parmread @@ -79,6 +80,8 @@ C 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 @@ -101,7 +104,7 @@ 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) @@ -139,7 +142,29 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword 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" + 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) + endif 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) @@ -161,6 +186,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword 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) @@ -217,7 +243,50 @@ cfmc modecalc=10 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 (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 @@ -344,8 +413,8 @@ C 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 @@ -539,6 +608,7 @@ C include 'COMMON.BOUNDS' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.SHIELD' character*4 sequence(maxres) integer rescode double precision x(maxvar) @@ -580,6 +650,9 @@ C Read weights of the subsequent energy terms. 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) @@ -718,7 +791,7 @@ C 12/1/95 Added weight for the multi-body term WCORR 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, @@ -856,11 +929,12 @@ C & write (iout,*) 'FTORS',ftors 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 - do i=1,nres - thetabound(1,i)=0 - thetabound(2,i)=pi - enddo +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 @@ -887,17 +961,17 @@ C E=k*x**4 & for_thet_constr(i) enddo endif - do i=1,nthet_constr + 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 - do i=1,ntheta_constr - ii = itheta_constr(i) - thetabound(1,ii) = phi0(i)-drange(i) - thetabound(2,ii) = phi0(i)+drange(i) - enddo +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 @@ -975,7 +1049,7 @@ c---------------------- 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 @@ -995,6 +1069,7 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) 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 @@ -1047,6 +1122,7 @@ C initial geometry. return else call read_angles(inp,*36) + call chainbuild_extconf endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' @@ -1071,6 +1147,7 @@ C initial geometry. 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.' @@ -1107,18 +1184,7 @@ C 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 @@ -1949,12 +2015,18 @@ C Get parameter filenames and open the parameter files. 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') @@ -1977,6 +2049,10 @@ c print *,"Processor",myrank," opened file IROTAM" 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') @@ -1989,6 +2065,8 @@ c print *,"Processor",myrank," opened file IFOURIER" 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) @@ -2006,6 +2084,10 @@ C Get parameter filenames and open the parameter files. 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) @@ -2014,6 +2096,8 @@ C Get parameter filenames and open the parameter files. 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) @@ -2030,6 +2114,10 @@ C Get parameter filenames and open the parameter files. 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 @@ -2044,11 +2132,15 @@ C Get parameter filenames and open the parameter files. 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 @@ -2118,7 +2210,7 @@ c print *,"Processor",myrank," fg_rank",fg_rank 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 @@ -2248,6 +2340,7 @@ c------------------------------------------------------------------------------- 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 @@ -2303,6 +2396,36 @@ c------------------------------------------------------------------------------- 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)