X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=e40717f716adead8d16ec6de4d313f0d555dd518;hb=36865e3203724663beeb1f60abdbafbec07b1019;hp=1cda5eac13362dfce355bd88dea59760e7c1eb10;hpb=c65d296f3cdc480d3737eae335a510332b4f0dc8;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 1cda5ea..e40717f 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 @@ -96,8 +99,12 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file 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) @@ -135,7 +142,24 @@ C Set up the time limit (caution! The time must be input in minutes!) 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,'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) @@ -157,6 +181,7 @@ C Set up the time limit (caution! The time must be input in minutes!) 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) @@ -213,7 +238,60 @@ 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 (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 @@ -340,8 +418,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 @@ -535,6 +613,7 @@ C include 'COMMON.BOUNDS' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.SHIELD' character*4 sequence(maxres) integer rescode double precision x(maxvar) @@ -576,6 +655,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) @@ -714,7 +796,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, @@ -851,6 +933,55 @@ C & write (iout,*) 'FTORS',ftors 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 @@ -923,7 +1054,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 @@ -943,6 +1074,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 @@ -995,6 +1127,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.' @@ -1019,6 +1152,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.' @@ -1055,18 +1189,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 @@ -1897,12 +2020,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') @@ -1925,6 +2054,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') @@ -1937,6 +2070,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) @@ -1954,6 +2089,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) @@ -1962,6 +2101,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) @@ -1978,6 +2119,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 @@ -1992,11 +2137,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 @@ -2066,7 +2215,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 @@ -2196,6 +2345,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 @@ -2251,6 +2401,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)