subroutine readrtns implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' include 'COMMON.SPLITELE' integer i,j logical file_exist C Read job setup parameters call read_control C Read force-field parameters except weights call parmread C Read control parameters for energy minimzation if required if (minim) call read_minim C Read MCM control parameters if required if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread C Read MD control parameters if reqjuired if (modecalc.eq.12) call read_MDpar C Read MREMD control parameters if required if (modecalc.eq.14) then call read_MDpar call read_REMDpar endif C Read MUCA control parameters if required if (lmuca) call read_muca C Read CSA control parameters if required (from fort.40 if exists C otherwise from general input file) if (modecalc.eq.8) then inquire (file="fort.40",exist=file_exist) if (.not.file_exist) call csaread endif cfmc if (modecalc.eq.10) call mcmfread c Read energy-term weights and disulfide parameters call weightread C Read molecule information, molecule geometry, energy-term weights, and C restraints if requested call molread C Print restraint information #ifdef MPI if (.not. out1file .or. me.eq.king) then #endif if (nhpb.gt.nss) then write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)') & "The following",nhpb-nss, & " distance restraints have been imposed:", & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V", & " score"," type" do i=nss+1,nhpb write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i), & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i), & irestr_type(i) enddo endif if (npeak.gt.0) then write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)') & "The following",npeak, & " NMR peak restraints have been imposed:", & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V", & " score"," type"," ipeak" do i=1,npeak do j=ipeak(1,i),ipeak(2,i) write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j), & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j), & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j) enddo enddo write (iout,*) "The ipeak array" do i=1,npeak write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i) enddo endif #ifdef MPI endif #endif c print *,"Processor",myrank," leaves READRTNS" return end C------------------------------------------------------------------------------- subroutine read_control C C Read contorl data C implicit none include 'DIMENSIONS' #ifdef MP include 'mpif.h' logical OKRandom, prng_restart double precision r1 #endif include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.THREAD' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.SAXS' include 'COMMON.MCM' include 'COMMON.MAP' include 'COMMON.HEADER' include 'COMMON.CSA' include 'COMMON.CHAIN' include 'COMMON.MUCA' include 'COMMON.MD' include 'COMMON.FFIELD' include 'COMMON.INTERACT' include 'COMMON.SETUP' include 'COMMON.SPLITELE' include 'COMMON.SHIELD' include 'COMMON.GEO' integer i integer KDIAG,ICORFL,IXDR COMMON /MACHSW/ KDIAG,ICORFL,IXDR character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ character*80 ucase character*320 controlcard double precision seed nglob_csa=0 eglob_csa=1d99 nmin_csa=0 read (INP,'(a)') titel call card_concat(controlcard) c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file call reada(controlcard,'SEED',seed,0.0D0) call random_init(seed) C Set up the time limit (caution! The time must be input in minutes!) read_cart=index(controlcard,'READ_CART').gt.0 out_cart=index(controlcard,'OUT_CART').gt.0 out_int=index(controlcard,'OUT_INT').gt.0 gmatout=index(controlcard,'GMATOUT').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,'NSAXS',nsaxs,0) call readi(controlcard,'SAXS_MODE',saxs_mode,0) call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0) call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0) write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE", & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff 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 call readi(controlcard,'SYM',symetr,1) 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 c call reada(controlcard,'RMSDBC',rmsdbc,3.0D0) c call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0) c call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0) c call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0) c call reada(controlcard,'DRMS',drms,0.1D0) c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then c write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc c write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 c write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max c write (iout,'(a,f10.1)')'DRMS = ',drms cc write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm c write (iout,'(a,f10.1)') 'Time limit (min):',timlim c endif call readi(controlcard,'NZ_START',nz_start,0) call readi(controlcard,'NZ_END',nz_end,0) c call readi(controlcard,'IZ_SC',iz_sc,0) timlim=60.0D0*timlim safety = 60.0d0*safety modecalc=0 call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100) call reada(controlcard,"T_BATH",t_bath,300.0d0) minim=(index(controlcard,'MINIMIZE').gt.0) dccart=(index(controlcard,'CART').gt.0) overlapsc=(index(controlcard,'OVERLAP').gt.0) overlapsc=.not.overlapsc searchsc=(index(controlcard,'SEARCHSC').gt.0 .and. & index(controlcard,'NOSEARCHSC').eq.0) sideadd=(index(controlcard,'SIDEADD').gt.0) energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) mremd_dec=(index(controlcard,'MREMD_DEC').gt.0) outpdb=(index(controlcard,'PDBOUT').gt.0) outmol2=(index(controlcard,'MOL2OUT').gt.0) pdbref=(index(controlcard,'PDBREF').gt.0) 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')) c print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'TUBEMOD',tubelog,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) call readi(controlcard,"RESCALE_MODE",rescale_mode,2) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) & write (iout,*) "RESCALE_MODE",rescale_mode split_ene=index(controlcard,'SPLIT_ENE').gt.0 if (index(controlcard,'REGULAR').gt.0.0D0) then call reada(controlcard,'WEIDIS',weidis,0.1D0) modecalc=1 refstr=.true. endif if (refstr) then if (index(controlcard,"CASC").gt.0) then iz_sc=1 c else if (index(controlcard,"CAONLY").gt.0) then c iz_sc=0 else if (index(controlcard,"SCONLY").gt.0) then iz_sc=2 else iz_sc=0 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when", c & " specifying REFSTR, PDBREF or REGULAR." c stop endif endif if (index(controlcard,'CHECKGRAD').gt.0) then modecalc=5 if (index(controlcard,' CART').gt.0) then icheckgrad=1 elseif (index(controlcard,'CARINT').gt.0) then icheckgrad=2 else icheckgrad=3 endif call reada(controlcard,'DELTA',aincr,1.0d-7) c write (iout,*) "icheckgrad",icheckgrad elseif (index(controlcard,'THREAD').gt.0) then modecalc=2 call readi(controlcard,'THREAD',nthread,0) if (nthread.gt.0) then call reada(controlcard,'WEIDIS',weidis,0.1D0) else if (fg_rank.eq.0) & write (iout,'(a)')'A number has to follow the THREAD keyword.' stop 'Error termination in Read_Control.' endif else if (index(controlcard,'MCMA').gt.0) then modecalc=3 else if (index(controlcard,'MCEE').gt.0) then modecalc=6 else if (index(controlcard,'MULTCONF').gt.0) then modecalc=4 else if (index(controlcard,'MAP').gt.0) then modecalc=7 call readi(controlcard,'MAP',nmap,0) else if (index(controlcard,'CSA').gt.0) then modecalc=8 crc else if (index(controlcard,'ZSCORE').gt.0) then crc crc ZSCORE is rm from UNRES, modecalc=9 is available crc crc modecalc=9 cfcm else if (index(controlcard,'MCMF').gt.0) then cfmc modecalc=10 else if (index(controlcard,'SOFTREG').gt.0) then modecalc=11 else if (index(controlcard,'CHECK_BOND').gt.0) then modecalc=-1 else if (index(controlcard,'TEST').gt.0) then modecalc=-2 else if (index(controlcard,'MD').gt.0) then modecalc=12 else if (index(controlcard,'RE ').gt.0) then modecalc=14 endif lmuca=index(controlcard,'MUCA').gt.0 call readi(controlcard,'MUCADYN',mucadyn,0) call readi(controlcard,'MUCASMOOTH',muca_smooth,0) if (lmuca .and. (me.eq.king .or. .not.out1file )) & then write (iout,*) 'MUCADYN=',mucadyn write (iout,*) 'MUCASMOOTH=',muca_smooth endif iscode=index(controlcard,'ONE_LETTER') indphi=index(controlcard,'PHI') indback=index(controlcard,'BACK') iranconf=index(controlcard,'RAND_CONF') 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) write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize c Cutoff range for interactions call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0) call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) write (iout,*) "Cutoff on interactions",r_cut_int write (iout,*) & "Cutoff in switching short and long range interactions in RESPA", & r_cut_respa write (iout,*) "lambda in switch function",rlamb 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 (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax if(me.eq.king.or..not.out1file) & write (iout,'(2a)') diagmeth(kdiag), & ' routine used to diagonalize matrices.' return end c-------------------------------------------------------------------------- subroutine read_REMDpar C C Read REMD settings C implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.MD' #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #endif include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.GEO' include 'COMMON.REMD' include 'COMMON.CONTROL' include 'COMMON.SETUP' character*80 ucase character*320 controlcard character*3200 controlcard1 integer iremd_m_total,i if(me.eq.king.or..not.out1file) & write (iout,*) "REMD setup" call card_concat(controlcard) call readi(controlcard,"NREP",nrep,3) call readi(controlcard,"NSTEX",nstex,1000) call reada(controlcard,"RETMIN",retmin,10.0d0) call reada(controlcard,"RETMAX",retmax,1000.0d0) mremdsync=(index(controlcard,'SYNC').gt.0) call readi(controlcard,"NSYN",i_sync_step,100) restart1file=(index(controlcard,'REST1FILE').gt.0) traj1file=(index(controlcard,'TRAJ1FILE').gt.0) call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1) if(max_cache_traj_use.gt.max_cache_traj) & max_cache_traj_use=max_cache_traj if(me.eq.king.or..not.out1file) then cd if (traj1file) then crc caching is in testing - NTWX is not ignored cd write (iout,*) "NTWX value is ignored" cd write (iout,*) " trajectory is stored to one file by master" cd write (iout,*) " before exchange at NSTEX intervals" cd endif write (iout,*) "NREP= ",nrep write (iout,*) "NSTEX= ",nstex write (iout,*) "SYNC= ",mremdsync write (iout,*) "NSYN= ",i_sync_step write (iout,*) "TRAJCACHE= ",max_cache_traj_use endif remd_tlist=.false. if (index(controlcard,'TLIST').gt.0) then remd_tlist=.true. call card_concat(controlcard1) read(controlcard1,*) (remd_t(i),i=1,nrep) if(me.eq.king.or..not.out1file) & write (iout,*)'tlist',(remd_t(i),i=1,nrep) endif remd_mlist=.false. if (index(controlcard,'MLIST').gt.0) then remd_mlist=.true. call card_concat(controlcard1) read(controlcard1,*) (remd_m(i),i=1,nrep) if(me.eq.king.or..not.out1file) then write (iout,*)'mlist',(remd_m(i),i=1,nrep) iremd_m_total=0 do i=1,nrep iremd_m_total=iremd_m_total+remd_m(i) enddo write (iout,*) 'Total number of replicas ',iremd_m_total endif endif if (adaptive) then write (iout,*) & "Adaptive (PMF-biased) umbrella sampling will be run" call PMFread endif if(me.eq.king.or..not.out1file) & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup " return end c-------------------------------------------------------------------------- subroutine read_MDpar C C Read MD settings C implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.MD' include 'COMMON.QRESTR' #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #endif include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.GEO' include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.SPLITELE' include 'COMMON.FFIELD' character*80 ucase character*320 controlcard integer i double precision eta call card_concat(controlcard) call readi(controlcard,"NSTEP",n_timestep,1000000) call readi(controlcard,"NTWE",ntwe,100) call readi(controlcard,"NTWX",ntwx,1000) call reada(controlcard,"DT",d_time,1.0d-1) call reada(controlcard,"DVMAX",dvmax,2.0d1) call reada(controlcard,"DAMAX",damax,1.0d1) call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1) call readi(controlcard,"LANG",lang,0) RESPA = index(controlcard,"RESPA") .gt. 0 call readi(controlcard,"NTIME_SPLIT",ntime_split,1) ntime_split0=ntime_split call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64) ntime_split0=ntime_split 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 scale_umb = index(controlcard,"SCALE_UMB").gt.0 adaptive = index(controlcard,"ADAPTIVE").gt.0 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) call reada(controlcard,"EQ_TIME",eq_time,1.0d+4) call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000) if (count_reset_moment.eq.0) count_reset_moment=1000000000 call readi(controlcard,"RESET_VEL",count_reset_vel,1000) reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0 if (count_reset_vel.eq.0) count_reset_vel=1000000000 large = index(controlcard,"LARGE").gt.0 print_compon = index(controlcard,"PRINT_COMPON").gt.0 rattle = index(controlcard,"RATTLE").gt.0 preminim = index(controlcard,"PREMINIM").gt.0 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then preminim=.true. dccart=.true. endif write (iout,*) "PREMINIM ",preminim if (preminim) then if (index(controlcard,'CART').gt.0) dccart=.true. write (iout,*) "dccart ",dccart write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR') if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim endif c if performing umbrella sampling, fragments constrained are read from the fragment file nset=0 if(usampl) then write (iout,*) "Umbrella sampling will be run" if (scale_umb.and.adaptive) then write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive" write (iout,*) "Select one of those and re-run the job." stop endif if (scale_umb) write (iout,*) &"Umbrella-restraint force constants will be scaled by temperature" call read_fragments endif if(me.eq.king.or..not.out1file) then write (iout,*) write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run " write (iout,*) write (iout,'(a)') "The units are:" write (iout,'(a)') "positions: angstrom, time: 48.9 fs" write (iout,'(2a)') "velocity: angstrom/(48.9 fs),", & " acceleration: angstrom/(48.9 fs)**2" write (iout,'(a)') "energy: kcal/mol, temperature: K" write (iout,*) write (iout,'(a60,i10)') "Number of time steps:",n_timestep write (iout,'(a60,f10.5,a)') & "Initial time step of numerical integration:",d_time, & " natural units" write (iout,'(60x,f10.5,a)') d_time*48.9," fs" write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int, & " A" write(iout,'(a60,i5)')"Frequency of updating interaction list", & imatupdate if (RESPA) then write (iout,'(2a,i4,a)') & "A-MTS algorithm used; initial time step for fast-varying", & " short-range forces split into",ntime_split," steps." write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff", & r_cut_respa," lambda",rlamb endif write (iout,'(2a,f10.5)') & "Maximum acceleration threshold to reduce the time step", & "/increase split number:",damax write (iout,'(2a,f10.5)') & "Maximum predicted energy drift to reduce the timestep", & "/increase split number:",edriftmax write (iout,'(a60,f10.5)') & "Maximum velocity threshold to reduce velocities:",dvmax write (iout,'(a60,i10)') "Frequency of property output:",ntwe write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx if (rattle) write (iout,'(a60)') & "Rattle algorithm used to constrain the virtual bonds" endif reset_fricmat=1000 if (lang.gt.0) then call reada(controlcard,"ETAWAT",etawat,0.8904d0) call reada(controlcard,"RWAT",rwat,1.4d0) call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2) surfarea=index(controlcard,"SURFAREA").gt.0 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000) if(me.eq.king.or..not.out1file)then write (iout,'(/a,$)') "Langevin dynamics calculation" if (lang.eq.1) then write (iout,'(a/)') & " with direct integration of Langevin equations" else if (lang.eq.2) then write (iout,'(a/)') " with TINKER stochasic MD integrator" else if (lang.eq.3) then write (iout,'(a/)') " with Ciccotti's stochasic MD integrator" else if (lang.eq.4) then write (iout,'(a/)') " in overdamped mode" else write (iout,'(//a,i5)') & "=========== ERROR: Unknown Langevin dynamics mode:",lang stop endif write (iout,'(a60,f10.5)') "Temperature:",t_bath write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat write (iout,'(a60,f10.5)') & "Scaling factor of the friction forces:",scal_fric if (surfarea) write (iout,'(2a,i10,a)') & "Friction coefficients will be scaled by solvent-accessible", & " surface area every",reset_fricmat," steps." endif c Calculate friction coefficients and bounds of stochastic forces eta=6*pi*cPoise*etawat if(me.eq.king.or..not.out1file) & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:" & ,eta gamp=scal_fric*(pstok+rwat)*eta stdfp=dsqrt(2*Rb*t_bath/d_time) do i=1,ntyp gamsc(i)=scal_fric*(restok(i)+rwat)*eta stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) enddo if(me.eq.king.or..not.out1file)then write (iout,'(/2a/)') & "Radii of site types and friction coefficients and std's of", & " stochastic forces of fully exposed sites" write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp) do i=1,ntyp write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i), & gamsc(i),stdfsc(i)*dsqrt(gamsc(i)) enddo endif else if (tbf) then if(me.eq.king.or..not.out1file)then write (iout,'(a)') "Berendsen bath calculation" write (iout,'(a60,f10.5)') "Temperature:",t_bath write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath if (reset_moment) & write (iout,'(a,i10,a)') "Momenta will be reset at zero every", & count_reset_moment," steps" if (reset_vel) & write (iout,'(a,i10,a)') & "Velocities will be reset at random every",count_reset_vel, & " steps" endif else if(me.eq.king.or..not.out1file) & write (iout,'(a31)') "Microcanonical mode calculation" endif if(me.eq.king.or..not.out1file)then if (rest) write (iout,'(/a/)') "===== Calculation restarted ====" if (usampl) then write(iout,*) "MD running with constraints." write(iout,*) "Equilibration time ", eq_time, " mtus." write(iout,*) "Constraining ", nfrag," fragments." write(iout,*) "Length of each fragment, weight and q0:" do iset=1,nset write (iout,*) "Set of restraints #",iset do i=1,nfrag write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset), & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset) enddo write(iout,*) "constraints between ", npair, "fragments." write(iout,*) "constraint pairs, weights and q0:" do i=1,npair write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset), & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset) enddo write(iout,*) "angle constraints within ", nfrag_back, & "backbone fragments." if (loc_qlike) then write(iout,*) "fragment, weights, q0:" do i=1,nfrag_back write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset), & ifrag_back(2,i,iset), & wfrag_back(1,i,iset),qin_back(1,i,iset), & wfrag_back(2,i,iset),qin_back(2,i,iset), & wfrag_back(3,i,iset),qin_back(3,i,iset) enddo else write(iout,*) "fragment, weights:" do i=1,nfrag_back write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset), & ifrag_back(2,i,iset),wfrag_back(1,i,iset), & wfrag_back(2,i,iset),wfrag_back(3,i,iset) enddo endif enddo iset=mod(kolor,nset)+1 endif endif if(me.eq.king.or..not.out1file) & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup " return end c------------------------------------------------------------------------------ subroutine molread C C Read molecular data. C implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer error_msg,ierror,ierr,ierrcode #endif include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.SAXS' include 'COMMON.DBASE' include 'COMMON.THREAD' include 'COMMON.CONTACTS' include 'COMMON.TORCNSTR' include 'COMMON.TIME1' include 'COMMON.BOUNDS' include 'COMMON.MD' include 'COMMON.SETUP' include 'COMMON.SHIELD' character*4 sequence(maxres) integer rescode double precision x(maxvar) character*256 pdbfile character*400 weightcard character*80 weightcard_t,ucase integer itype_pdb(maxres) common /pizda/ itype_pdb logical seq_comp,fail double precision energia(0:n_ene) double precision secprob(3,maxdih_constr) double precision co double precision phihel,phibet,sigmahel,sigmabet integer iti,nsi,maxsi integer ilen external ilen integer iperm,tperm integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2 double precision sumv C C Read PDB structure if applicable C if (indpdb.gt.0 .or. pdbref) then read(inp,'(a)') pdbfile if(me.eq.king.or..not.out1file) & write (iout,'(2a)') 'PDB data will be read from file ', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a)') 'Error opening PDB file.' stop 34 continue call readpdb do i=1,2*nres do j=1,3 crefjlee(j,i)=c(j,i) enddo enddo #ifdef DEBUG do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo #endif if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup, & ' nstart_sup=',nstart_sup do i=1,nres itype_pdb(i)=itype(i) enddo close (ipdbin) nnt=nstart_sup nct=nstart_sup+nsup-1 call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then if(me.eq.king.or..not.out1file) & write(iout,*)'Adding sidechains' maxsi=1000 do i=2,nres-1 iti=itype(i) if (iti.ne.10 .and. itype(i).ne.ntyp1) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) call gen_side(iti,theta(i+1),alph(i),omeg(i),fail) nsi=nsi+1 enddo if(fail) write(iout,*)'Adding sidechain failed for res ', & i,' after ',nsi,' trials' endif enddo endif endif if (indpdb.eq.0) then C Read sequence if not taken from the pdb file. read (inp,*) nres c print *,'nres=',nres if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) else read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) endif C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo endif c print *,nres c print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR else if (iabs(itype(i+1)).ne.20) then #else else if (iabs(itype(i)).ne.20) then #endif itel(i)=1 else itel(i)=2 endif enddo if(me.eq.king.or..not.out1file)then write (iout,*) "ITEL" do i=1,nres-1 write (iout,*) i,itype(i),itel(i) enddo print *,'Call Read_Bridge.' endif nnt=1 nct=nres cd print *,'NNT=',NNT,' NCT=',NCT call seq2chains(nres,itype,nchain,chain_length,chain_border, & ireschain) chain_border1(1,1)=1 chain_border1(2,1)=chain_border(2,1)+1 do i=2,nchain-1 chain_border1(1,i)=chain_border(1,i)-1 chain_border1(2,i)=chain_border(2,i)+1 enddo chain_border1(1,nchain)=chain_border(1,nchain)-1 chain_border1(2,nchain)=nres write(iout,*) "nres",nres," nchain",nchain do i=1,nchain write(iout,*)"chain",i,chain_length(i),chain_border(1,i), & chain_border(2,i),chain_border1(1,i),chain_border1(2,i) enddo call chain_symmetry(nchain,nres,itype,chain_border, & chain_length,npermchain,tabpermchain) c do i=1,nres c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain), c & ii=1,npermchain) c enddo write(iout,*) "residue permutations" do i=1,nres write(iout,*) i,(iperm(i,ii),ii=1,npermchain) enddo call flush(iout) if (itype(1).eq.ntyp1) nnt=2 if (itype(nres).eq.ntyp1) nct=nct-1 write (iout,*) "nnt",nnt," nct",nct call flush(iout) #ifdef DFA if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and. & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then call init_dfa_vars print*, 'init_dfa_vars finished!' call read_dfa_info print*, 'read_dfa_info finished!' endif #endif if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup nstart_seq=nnt if (nsup.le.(nct-nnt+1)) then do i=0,nct-nnt+1-nsup if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then nstart_seq=nnt+i goto 111 endif enddo write (iout,'(a)') & 'Error - sequences to be superposed do not match.' stop else do i=0,nsup-(nct-nnt+1) if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) & then nstart_sup=nstart_sup+i nsup=nct-nnt+1 goto 111 endif enddo write (iout,'(a)') & 'Error - sequences to be superposed do not match.' endif 111 continue if (nsup.eq.0) nsup=nct-nnt if (nstart_sup.eq.0) nstart_sup=nnt if (nstart_seq.eq.0) nstart_seq=nnt if(me.eq.king.or..not.out1file) & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup, & ' nstart_seq=',nstart_seq endif call read_bridge C 8/13/98 Set limits to generating the dihedral angles do i=1,nres phibound(1,i)=-pi phibound(2,i)=pi enddo read (inp,*) ndih_constr write (iout,*) "ndih_constr",ndih_constr if (ndih_constr.gt.0) then raw_psipred=.false. C read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), & i=1,ndih_constr) #ifdef MPI if(me.eq.king.or..not.out1file)then #endif write (iout,*) & 'There are',ndih_constr,' restraints on gamma angles.' do i=1,ndih_constr write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), & ftors(i) enddo #ifdef MPI endif #endif do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo C if(me.eq.king.or..not.out1file) C & write (iout,*) 'FTORS',ftors do i=1,ndih_constr ii = idih_constr(i) phibound(1,ii) = phi0(i)-drange(i) phibound(2,ii) = phi0(i)+drange(i) enddo #ifdef MPI if (me.eq.king .or. .not.out1file) then #endif write (iout,'(a)') 'Boundaries in gamma angle sampling:' do i=1,nres write (iout,'(a3,i5,2f10.1)') & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg enddo #ifdef MP endif #endif else if (ndih_constr.lt.0) then raw_psipred=.true. call card_concat(weightcard) 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).ne.ntyp1 .and. itype(i-2).ne.ntyp1 & .and. itype(i-1).ne.ntyp1 .and. itype(i).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 #ifdef MPI if(me.eq.king.or..not.out1file)then #endif write (iout,*) & 'There are',ndih_constr, & ' bimodal restraints on gamma angles.' do i=1,ndih_constr write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i, & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2, & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1, & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg, & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg, & (vpsipred(j,i),j=1,3) enddo #ifdef MPI endif #endif c Raw psipred input 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 c--- Zscore rms ------- if (nz_start.eq.0) nz_start=nnt if (nz_end.eq.0 .and. nsup.gt.0) then nz_end=nnt+nsup-1 else if (nz_end.eq.0) then nz_end=nct endif if(me.eq.king.or..not.out1file)then write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end write (iout,*) 'IZ_SC=',iz_sc endif c---------------------- call init_int_table if (refstr) then if (.not.pdbref) then call read_angles(inp,*38) call bond_regular goto 39 38 write (iout,'(a)') 'Error reading reference structure.' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,IERROR) stop 'Error reading reference structure' #endif 39 call chainbuild_extconf call setup_var czscore call geom_to_var(nvar,coord_exp_zs(1,1)) nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 do i=1,2*nres do j=1,3 cref(j,i)=c(j,i) enddo enddo call contact(.true.,ncont_ref,icont_ref,co) endif endif c print *, "A TU" write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup 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 call NMRpeak_partition if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co if (pdbref) then if(me.eq.king.or..not.out1file) & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup do i=1,ncont_ref do j=1,2 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup enddo if(me.eq.king.or..not.out1file) & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ', & icont_ref(1,i),' ', & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) enddo endif write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr if (constr_homology.gt.0) then call read_constr_homology if (indpdb.gt.0 .or. pdbref) then do i=1,2*nres do j=1,3 c(j,i)=crefjlee(j,i) cref(j,i)=crefjlee(j,i) enddo enddo endif #ifdef DEBUG write (iout,*) "sc_loc_geom: Array C" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3), & (c(j,i+nres),j=1,3) enddo write (iout,*) "Array Cref" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3), & (cref(j,i+nres),j=1,3) enddo #endif call int_from_cart1(.false.) 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 enddo else homol_nset=0 endif C endif if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then C If input structure hasn't been supplied from the PDB file read or generate C initial geometry. if (iranconf.eq.0 .and. .not. extconf) then if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) & write (iout,'(a)') 'Initial geometry will be read in.' if (read_cart) then read(inp,'(8f10.5)',end=36,err=36) & ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "Exit READ_CART" c write (iout,'(8f10.5)') c & ((c(l,k),l=1,3),k=1,nres), c & ((c(l,k+nres),l=1,3),k=nnt,nct) call cartprint do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1) enddo enddo do i=nnt,nct if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres) enddo else do j=1,3 dc(j,i+nres)=0.0d0 c dc_norm(j,i+nres)=0.0d0 enddo endif enddo call int_from_cart1(.true.) write (iout,*) "Finish INT_TO_CART" c write (iout,*) "DC" c do i=1,nres c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3), c & (dc(j,i+nres),j=1,3) c enddo c call cartprint c call setup_var c return else call read_angles(inp,*36) call bond_regular call chainbuild_extconf endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' #ifdef MPI call mpi_finalize( MPI_COMM_WORLD,IERR ) #endif stop 'Error reading angle file.' 37 continue else if (extconf) then if (me.eq.king.or..not.out1file .and. fg_rank.eq.0) & write (iout,'(a)') 'Extended chain initial geometry.' do i=3,nres theta(i)=90d0*deg2rad enddo do i=4,nres phi(i)=180d0*deg2rad enddo do i=2,nres-1 alph(i)=110d0*deg2rad enddo do i=2,nres-1 omeg(i)=-120d0*deg2rad if (itype(i).le.0) omeg(i)=-omeg(i) enddo call bond_regular call chainbuild_extconf else if(me.eq.king.or..not.out1file) & write (iout,'(a)') 'Random-generated initial geometry.' call bond_regular #ifdef MPI if (me.eq.king .or. fg_rank.eq.0 .and. ( & modecalc.eq.12 .or. modecalc.eq.14) ) then #endif 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 (*,*) 'Processor:',me, & ' Failed to generate random conformation', & ' itrial=',itrial call intout #ifdef AIX call flush_(iout) #else call flush(iout) #endif enddo write (iout,'(a,i3,a)') 'Processor:',me, & ' error in generating random conformation.' write (*,'(a,i3,a)') 'Processor:',me, & ' error in generating random conformation.' call flush(iout) #ifdef MPI call MPI_Abort(mpi_comm_world,error_msg,ierrcode) 40 continue endif #else write (*,'(a)') & ' error in generating random conformation.' stop 40 continue #endif endif elseif (modecalc.eq.4) then read (inp,'(a)') intinname open (intin,file=intinname,status='old',err=333) if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) & write (iout,'(a)') 'intinname',intinname write (*,'(a)') 'Processor',myrank,' intinname',intinname goto 334 333 write (iout,'(2a)') 'Error opening angle file ',intinname #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,IERR) #endif stop 'Error opening angle file.' 334 continue endif C Generate distance constraints, if the PDB structure is to be regularized. if (nthread.gt.0) then call read_threadbase endif call setup_var if (me.eq.king .or. .not. out1file) & call intout if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) if (dyn_ss) then write(iout,*)"Running with dynamic disulfide-bond formation" else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres it1=itype(i1) it2=itype(i2) write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), & ebr,forcon(i) enddo write (iout,'(a)') endif endif if (ns.gt.0.and.dyn_ss) then do i=nss+1,nhpb ihpb(i-nss)=ihpb(i) jhpb(i-nss)=jhpb(i) forcon(i-nss)=forcon(i) dhpb(i-nss)=dhpb(i) enddo nhpb=nhpb-nss nss=0 call hpb_partition do i=1,ns dyn_ss_mask(iss(i))=.true. enddo endif c --- test c call cartprint c write (iout,*) "DC" c do i=1,nres c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3), c & (dc(j,i+nres),j=1,3) c enddo if (i2ndstr.gt.0) call secstrp2dihc c call geom_to_var(nvar,x) c call etotal(energia(0)) c call enerprint(energia(0)) c call briefout(0,etot) c stop cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT cd write (iout,'(a)') 'Variable list:' cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar) #ifdef MPI if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') & 'Processor',myrank,': end reading molecular data.' #endif c print *,"A TU?" return end c-------------------------------------------------------------------------- logical function seq_comp(itypea,itypeb,length) implicit none integer length,itypea(length),itypeb(length) integer i do i=1,length if (itypea(i).ne.itypeb(i)) then seq_comp=.false. return endif enddo seq_comp=.true. return end c----------------------------------------------------------------------------- subroutine read_bridge C Read information about disulfide bridges. implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierror #endif include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.DBASE' include 'COMMON.THREAD' include 'COMMON.TIME1' include 'COMMON.SETUP' integer i,j C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) print *,'ns=',ns if(me.eq.king.or..not.out1file) & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns) C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ', & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ', & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) stop #endif endif enddo C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) if(fg_rank.eq.0) & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of C bridging residues. do i=1,nss do j=1,i-1 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then write (iout,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' write (*,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) stop #endif endif enddo do j=1,ns if (ihpb(i).eq.iss(j)) goto 10 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 10 continue do j=1,ns if (jhpb(i).eq.iss(j)) goto 20 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue c dhpb(i)=dbr c forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo endif endif return end c---------------------------------------------------------------------------- subroutine read_x(kanal,*) implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.LOCAL' include 'COMMON.INTERACT' integer i,j,k,l,kanal c Read coordinates from input c read(kanal,'(8f10.5)',end=10,err=10) & ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) do j=1,3 c(j,nres+1)=c(j,1) c(j,2*nres)=c(j,nres) enddo call int_from_cart1(.false.) 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=nnt,nct if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then 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 endif enddo return 10 return1 end c---------------------------------------------------------------------------- subroutine read_threadbase implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.DBASE' include 'COMMON.THREAD' include 'COMMON.TIME1' integer i,j,k double precision dist C Read pattern database for threading. read (icbase,*) nseq do i=1,nseq read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i), & nres_base(2,i),nres_base(3,i) read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1, & nres_base(1,i)) c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i), c & nres_base(2,i),nres_base(3,i) c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1, c & nres_base(1,i)) enddo close (icbase) if (weidis.eq.0.0D0) weidis=0.1D0 do i=nnt,nct do j=i+2,nct nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=weidis enddo enddo read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl) write (iout,'(a,i5)') 'nexcl: ',nexcl write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl) return end c------------------------------------------------------------------------------ subroutine setup_var implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.DBASE' include 'COMMON.THREAD' include 'COMMON.TIME1' integer i C Set up variable list. ntheta=nres-2 nphi=nres-3 nvar=ntheta+nphi nside=0 do i=2,nres-1 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then nside=nside+1 ialph(i,1)=nvar+nside ialph(nside,2)=i endif enddo if (indphi.gt.0) then nvar=nphi else if (indback.gt.0) then nvar=nphi+ntheta else nvar=nvar+2*nside endif return end c---------------------------------------------------------------------------- subroutine gen_dist_constr C Generate CA distance constraints. implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.DBASE' include 'COMMON.THREAD' include 'COMMON.TIME1' integer i,j,itype_pdb(maxres) common /pizda/ itype_pdb double precision dist character*2 iden cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct, cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq, cd & ' nsup',nsup 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)) do j=i+2,nstart_sup+nsup-1 nhpb=nhpb+1 ihpb(nhpb)=i+nstart_seq-nstart_sup jhpb(nhpb)=j+nstart_seq-nstart_sup forcon(nhpb)=weidis dhpb(nhpb)=dist(i,j) enddo enddo cd write (iout,'(a)') 'Distance constraints:' cd do i=nss+1,nhpb cd ii=ihpb(i) cd jj=jhpb(i) cd iden='CA' cd if (ii.gt.nres) then cd iden='SC' cd ii=ii-nres cd jj=jj-nres cd endif cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)') cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj, cd & dhpb(i),forcon(i) cd enddo return end c---------------------------------------------------------------------------- subroutine map_read implicit none include 'DIMENSIONS' include 'COMMON.MAP' include 'COMMON.IOUNITS' integer imap character*3 angid(4) /'THE','PHI','ALP','OME'/ character*80 mapcard,ucase do imap=1,nmap read (inp,'(a)') mapcard mapcard=ucase(mapcard) if (index(mapcard,'PHI').gt.0) then kang(imap)=1 else if (index(mapcard,'THE').gt.0) then kang(imap)=2 else if (index(mapcard,'ALP').gt.0) then kang(imap)=3 else if (index(mapcard,'OME').gt.0) then kang(imap)=4 else write(iout,'(a)')'Error - illegal variable spec in MAP card.' stop 'Error - illegal variable spec in MAP card.' endif call readi (mapcard,'RES1',res1(imap),0) call readi (mapcard,'RES2',res2(imap),0) if (res1(imap).eq.0) then res1(imap)=res2(imap) else if (res2(imap).eq.0) then res2(imap)=res1(imap) endif if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then write (iout,'(a)') & 'Error - illegal definition of variable group in MAP.' stop 'Error - illegal definition of variable group in MAP.' endif call reada(mapcard,'FROM',ang_from(imap),0.0D0) call reada(mapcard,'TO',ang_to(imap),0.0D0) call readi(mapcard,'NSTEP',nstep(imap),0) if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then write (iout,'(a)') & 'Illegal boundary and/or step size specification in MAP.' stop 'Illegal boundary and/or step size specification in MAP.' endif enddo ! imap return end c---------------------------------------------------------------------------- subroutine csaread implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.CSA' include 'COMMON.BANK' include 'COMMON.CONTROL' character*80 ucase character*620 mcmcard call card_concat(mcmcard) call readi(mcmcard,'NCONF',nconf,50) call readi(mcmcard,'NADD',nadd,0) call readi(mcmcard,'JSTART',jstart,1) call readi(mcmcard,'JEND',jend,1) call readi(mcmcard,'NSTMAX',nstmax,500000) call readi(mcmcard,'N0',n0,1) call readi(mcmcard,'N1',n1,6) call readi(mcmcard,'N2',n2,4) call readi(mcmcard,'N3',n3,0) call readi(mcmcard,'N4',n4,0) call readi(mcmcard,'N5',n5,0) call readi(mcmcard,'N6',n6,10) call readi(mcmcard,'N7',n7,0) call readi(mcmcard,'N8',n8,0) call readi(mcmcard,'N9',n9,0) call readi(mcmcard,'N14',n14,0) call readi(mcmcard,'N15',n15,0) call readi(mcmcard,'N16',n16,0) call readi(mcmcard,'N17',n17,0) call readi(mcmcard,'N18',n18,0) vdisulf=(index(mcmcard,'DYNSS').gt.0) call readi(mcmcard,'NDIFF',ndiff,2) call reada(mcmcard,'DIFFCUT',diffcut,0.0d0) call readi(mcmcard,'IS1',is1,1) call readi(mcmcard,'IS2',is2,8) call readi(mcmcard,'NRAN0',nran0,4) call readi(mcmcard,'NRAN1',nran1,2) call readi(mcmcard,'IRR',irr,1) call readi(mcmcard,'NSEED',nseed,20) call readi(mcmcard,'NTOTAL',ntotal,10000) call reada(mcmcard,'CUT1',cut1,2.0d0) call reada(mcmcard,'CUT2',cut2,5.0d0) call reada(mcmcard,'ESTOP',estop,-3000.0d0) call readi(mcmcard,'ICMAX',icmax,3) call readi(mcmcard,'IRESTART',irestart,0) c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0) ntbankm=0 c!bankt call reada(mcmcard,'DELE',dele,20.0d0) call reada(mcmcard,'DIFCUT',difcut,720.0d0) call readi(mcmcard,'IREF',iref,0) call reada(mcmcard,'RMSCUT',rmscut,4.0d0) call reada(mcmcard,'PNCCUT',pnccut,0.5d0) call readi(mcmcard,'NCONF_IN',nconf_in,0) call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0) write (iout,*) "NCONF_IN",nconf_in return end c---------------------------------------------------------------------------- subroutine mcmread implicit none include 'DIMENSIONS' include 'COMMON.MCM' include 'COMMON.MCE' include 'COMMON.IOUNITS' character*80 ucase character*320 mcmcard integer i call card_concat(mcmcard) call readi(mcmcard,'MAXACC',maxacc,100) call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000) call readi(mcmcard,'MAXTRIAL',maxtrial,100) call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000) call readi(mcmcard,'MAXREPM',maxrepm,200) call reada(mcmcard,'RANFRACT',RanFract,0.5D0) call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0) call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3) call reada(mcmcard,'E_UP',e_up,5.0D0) call reada(mcmcard,'DELTE',delte,0.1D0) call readi(mcmcard,'NSWEEP',nsweep,5) call readi(mcmcard,'NSTEPH',nsteph,0) call readi(mcmcard,'NSTEPC',nstepc,0) call reada(mcmcard,'TMIN',tmin,298.0D0) call reada(mcmcard,'TMAX',tmax,298.0D0) call readi(mcmcard,'NWINDOW',nwindow,0) call readi(mcmcard,'PRINT_MC',print_mc,0) print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0) print_int=(index(mcmcard,'NO_PRINT_INT').le.0) ent_read=(index(mcmcard,'ENT_READ').gt.0) call readi(mcmcard,'SAVE_FREQ',save_frequency,1000) call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000) call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000) call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000) call readi(mcmcard,'PRINT_FREQ',print_freq,1000) if (nwindow.gt.0) then read (inp,*) (winstart(i),winend(i),i=1,nwindow) do i=1,nwindow winlen(i)=winend(i)-winstart(i)+1 enddo endif if (tmax.lt.tmin) tmax=tmin if (tmax.eq.tmin) then nstepc=0 nsteph=0 endif if (nstepc.gt.0 .and. nsteph.gt.0) then tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0)) tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0)) endif C Probabilities of different move types sumpro_type(0)=0.0D0 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0) call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0) sumpro_type(2)=sumpro_type(1)+sumpro_type(2) call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0) sumpro_type(3)=sumpro_type(2)+sumpro_type(3) call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0) sumpro_type(4)=sumpro_type(3)+sumpro_type(4) do i=1,MaxMoveType print *,'i',i,' sumprotype',sumpro_type(i) sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType) print *,'i',i,' sumprotype',sumpro_type(i) enddo return end c---------------------------------------------------------------------------- subroutine read_minim implicit none include 'DIMENSIONS' include 'COMMON.MINIM' include 'COMMON.IOUNITS' character*80 ucase character*320 minimcard call card_concat(minimcard) call readi(minimcard,'MAXMIN',maxmin,2000) call readi(minimcard,'MAXFUN',maxfun,5000) call readi(minimcard,'MINMIN',minmin,maxmin) call readi(minimcard,'MINFUN',minfun,maxmin) call reada(minimcard,'TOLF',tolf,1.0D-2) call reada(minimcard,'RTOLF',rtolf,1.0D-4) print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1) print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1) print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1) write (iout,'(/80(1h*)/20x,a/80(1h*))') & 'Options in energy minimization:' write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') & 'MaxMin:',MaxMin,' MaxFun:',MaxFun, & 'MinMin:',MinMin,' MinFun:',MinFun, & ' TolF:',TolF,' RTolF:',RTolF return end c---------------------------------------------------------------------------- subroutine read_angles(kanal,*) implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' integer i,kanal c Read angles from input c read (kanal,*,err=10,end=10) (theta(i),i=3,nres) read (kanal,*,err=10,end=10) (phi(i),i=4,nres) read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1) read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1) do i=1,nres c 9/7/01 avoid 180 deg valence angle if (theta(i).gt.179.99d0) theta(i)=179.99d0 c theta(i)=deg2rad*theta(i) phi(i)=deg2rad*phi(i) alph(i)=deg2rad*alph(i) omeg(i)=deg2rad*omeg(i) enddo return 10 return1 end c---------------------------------------------------------------------------- subroutine reada(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch double precision wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,err=10,end=10) wartosc return 10 wartosc=default return end c---------------------------------------------------------------------------- subroutine readi(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch integer wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,err=10,end=10) wartosc return 10 wartosc=default return end c---------------------------------------------------------------------------- subroutine multreadi(rekord,lancuch,tablica,dim,default) implicit none integer dim,i integer tablica(dim),default character*(*) rekord,lancuch character*80 aux integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine multreada(rekord,lancuch,tablica,dim,default) implicit none integer dim,i double precision tablica(dim),default character*(*) rekord,lancuch character*80 aux integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine openunits implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierror character*16 form,nodename integer nodelen #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' include 'COMMON.MD' include 'COMMON.CONTROL' integer lenpre,lenpot,ilen,lentmp,npos external ilen character*3 out1file_text,ucase character*3 ll external ucase tmpdir="" out1file_text="YES" c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits" call getenv_loc("PREFIX",prefix) pref_orig = prefix call getenv_loc("POT",pot) call getenv_loc("DIRTMP",tmpdir) call getenv_loc("CURDIR",curdir) call getenv_loc("OUT1FILE",out1file_text) c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV" out1file_text=ucase(out1file_text) if (out1file_text(1:1).eq."Y") then out1file=.true. else out1file=fg_rank.gt.0 endif lenpre=ilen(prefix) lenpot=ilen(pot) lentmp=ilen(tmpdir) if (lentmp.gt.0) then write (*,'(80(1h!))') write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!" write (*,'(80(1h!))') write (*,*)"All output files will be on node /tmp directory." #ifdef MPI call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR ) if (me.eq.king) then write (*,*) "The master node is ",nodename else if (fg_rank.eq.0) then write (*,*) "I am the CG slave node ",nodename else write (*,*) "I am the FG slave node ",nodename endif #endif PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre) lenpre = lentmp+lenpre+1 endif entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr' C Get the names and open the input files #if defined(WINIFL) || defined(WINPGI) open(1,file=pref_orig(:ilen(pref_orig))// & '.inp',status='old',readonly,shared) open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown') C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown') C Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',readonly,shared) call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',readonly,shared) call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',readonly,shared) call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old',readonly,shared) call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,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 1" open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown') c print *,"Processor",myrank," opened file 9" C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown') C Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') c print *,"Processor",myrank," opened file IBOND" call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) open (ithep_pdb,file=thetname_pdb,status='old',action='read') #endif c print *,"Processor",myrank," opened file ITHEP" call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif c print *,"Processor",myrank," opened file IROTAM" call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old',action='read') c print *,"Processor",myrank," opened file ITORP" call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,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 ISCCOR" call getenv_loc('FOURIER',fouriername) open (ifourier,file=fouriername,status='old',action='read') c print *,"Processor",myrank," opened file IFOURIER" call getenv_loc('ELEPAR',elename) open (ielep,file=elename,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(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old') open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown') C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown') C Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old') call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old') #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) open (ithep_pdb,file=thetname_pdb,status='old') #endif call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old') #endif call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old') call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old') call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old') call getenv_loc('FOURIER',fouriername) open (ifourier,file=fouriername,status='old') call getenv_loc('ELEPAR',elename) 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 (9,file=prefix(:ilen(prefix))//'.intin',status='unknown') C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown') C Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',readonly) call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',readonly) call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',readonly) call getenv_loc('TORPAR',torname) open (itorp,file=torname,status='old',readonly) call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old',readonly) call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old',readonly) #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) open (ithep_pdb,file=thetname_pdb,status='old',action='read') #endif call getenv_loc('FOURIER',fouriername) open (ifourier,file=fouriername,status='old',readonly) call getenv_loc('ELEPAR',elename) 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 #ifdef TUBE call getenv_loc('TUBEPAR',tubename) #if defined(WINIFL) || defined(WINPGI) open (itube,file=tubename,status='old',readonly,shared) #elif (defined CRAY) || (defined AIX) open (itube,file=tubename,status='old',action='read') #elif (defined G77) open (itube,file=tubename,status='old') #else open (itube,file=tubename,status='old',readonly) #endif #endif #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file C Use -DOLDSCP to use hard-coded constants instead. C call getenv_loc('SCPPAR',scpname) #if defined(WINIFL) || defined(WINPGI) open (iscpp,file=scpname,status='old',readonly,shared) #elif (defined CRAY) || (defined AIX) open (iscpp,file=scpname,status='old',action='read') #elif (defined G77) open (iscpp,file=scpname,status='old') #else open (iscpp,file=scpname,status='old',readonly) #endif #endif call getenv_loc('PATTERN',patname) #if defined(WINIFL) || defined(WINPGI) open (icbase,file=patname,status='old',readonly,shared) #elif (defined CRAY) || (defined AIX) open (icbase,file=patname,status='old',action='read') #elif (defined G77) open (icbase,file=patname,status='old') #else open (icbase,file=patname,status='old',readonly) #endif #ifdef MPI C Open output file only for CG processes c print *,"Processor",myrank," fg_rank",fg_rank if (fg_rank.eq.0) then if (nodes.eq.1) then npos=3 else npos = dlog10(dfloat(nodes-1))+1 endif if (npos.lt.3) npos=3 write (liczba,'(i1)') npos form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) & //')' write (liczba,form) me outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// & liczba(:ilen(liczba)) intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) & //'.int' pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) & //'.pdb' mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.mol2' statname=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.stat' if (lentmp.gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) & //liczba(:ilen(liczba))//'.stat') rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) & //'.rst' if(usampl) then qname=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.const' endif endif #else outname=prefix(:lenpre)//'.out_'//pot(:lenpot) intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int' pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb' mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2' statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat' if (lentmp.gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) & //'.stat') rest2name=prefix(:ilen(prefix))//'.rst' if(usampl) then qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const' endif #endif #if defined(AIX) || defined(PGI) || defined(CRAY) if (me.eq.king .or. .not. out1file) then open(iout,file=outname,status='unknown') else open(iout,file="/dev/null",status="unknown") endif c#define DEBUG #ifdef DEBUG if (fg_rank.gt.0) then write (liczba,'(i3.3)') myrank/nfgtasks write (ll,'(bz,i3.3)') fg_rank open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll, & status='unknown') endif #endif c#undef DEBUG if(me.eq.king) then open(igeom,file=intname,status='unknown',position='append') open(ipdb,file=pdbname,status='unknown') open(imol2,file=mol2name,status='unknown') open(istat,file=statname,status='unknown',position='append') else c1out open(iout,file=outname,status='unknown') endif #else if (me.eq.king .or. .not.out1file) & open(iout,file=outname,status='unknown') #ifdef DEBUG if (fg_rank.gt.0) then write (liczba,'(i3.3)') myrank/nfgtasks write (ll,'(bz,i3.3)') fg_rank open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll, & status='unknown') endif #endif if(me.eq.king) then open(igeom,file=intname,status='unknown',access='append') open(ipdb,file=pdbname,status='unknown') open(imol2,file=mol2name,status='unknown') open(istat,file=statname,status='unknown',access='append') else c1out open(iout,file=outname,status='unknown') endif #endif csa_rbank=prefix(:lenpre)//'.CSA.rbank' csa_seed=prefix(:lenpre)//'.CSA.seed' csa_history=prefix(:lenpre)//'.CSA.history' csa_bank=prefix(:lenpre)//'.CSA.bank' csa_bank1=prefix(:lenpre)//'.CSA.bank1' csa_alpha=prefix(:lenpre)//'.CSA.alpha' csa_alpha1=prefix(:lenpre)//'.CSA.alpha1' c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt' csa_int=prefix(:lenpre)//'.int' csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized' csa_native_int=prefix(:lenpre)//'.CSA.native.int' csa_in=prefix(:lenpre)//'.CSA.in' c print *,"Processor",myrank,"fg_rank",fg_rank," opened files" C Write file names if (me.eq.king)then write (iout,'(80(1h-))') write (iout,'(30x,a)') "FILE ASSIGNMENT" write (iout,'(80(1h-))') write (iout,*) "Input file : ", & pref_orig(:ilen(pref_orig))//'.inp' write (iout,*) "Output file : ", & outname(:ilen(outname)) write (iout,*) write (iout,*) "Sidechain potential file : ", & sidename(:ilen(sidename)) #ifndef OLDSCP write (iout,*) "SCp potential file : ", & scpname(:ilen(scpname)) #endif write (iout,*) "Electrostatic potential file : ", & elename(:ilen(elename)) write (iout,*) "Cumulant coefficient file : ", & fouriername(:ilen(fouriername)) write (iout,*) "Torsional parameter file : ", & torname(:ilen(torname)) write (iout,*) "Double torsional parameter file : ", & tordname(:ilen(tordname)) write (iout,*) "SCCOR parameter file : ", & sccorname(:ilen(sccorname)) write (iout,*) "Bond & inertia constant file : ", & bondname(:ilen(bondname)) write (iout,*) "Bending-torsion parameter file : ", & thetname(:ilen(thetname)) write (iout,*) "Rotamer parameter file : ", & rotname(:ilen(rotname)) write (iout,*) "Thetpdb parameter file : ", & thetname_pdb(:ilen(thetname_pdb)) write (iout,*) "Threading database : ", & patname(:ilen(patname)) if (lentmp.ne.0) &write (iout,*)" DIRTMP : ", & tmpdir(:lentmp) write (iout,'(80(1h-))') endif return end c---------------------------------------------------------------------------- subroutine card_concat(card) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.IOUNITS' character*(*) card character*80 karta,ucase external ilen read (inp,'(a)') karta karta=ucase(karta) card=' ' do while (karta(80:80).eq.'&') card=card(:ilen(card)+1)//karta(:79) read (inp,'(a)') karta karta=ucase(karta) enddo card=card(:ilen(card)+1)//karta return end c------------------------------------------------------------------------------ subroutine readrst implicit none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.MD' include 'COMMON.QRESTR' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif integer i,j 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 do i=1,2*nres read(irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo if(usampl) then read (irest2,*) iset endif close(irest2) return end c------------------------------------------------------------------------------ subroutine read_fragments implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.MD' include 'COMMON.QRESTR' include 'COMMON.CONTROL' integer i read(inp,*) nset,nfrag,npair,nfrag_back loc_qlike=(nfrag_back.lt.0) nfrag_back=iabs(nfrag_back) if(me.eq.king.or..not.out1file) & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair, & " nfrag_back",nfrag_back," loc_qlike",loc_qlike do iset=1,nset read(inp,*) mset(iset) do i=1,nfrag read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset), & qinfrag(i,iset) if(me.eq.king.or..not.out1file) & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset), & ifrag(2,i,iset), qinfrag(i,iset) enddo do i=1,npair read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset), & qinpair(i,iset) if(me.eq.king.or..not.out1file) & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset), & ipair(2,i,iset), qinpair(i,iset) enddo if (loc_qlike) then do i=1,nfrag_back read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset), & wfrag_back(2,i,iset),qin_back(2,i,iset), & wfrag_back(3,i,iset),qin_back(3,i,iset), & ifrag_back(1,i,iset),ifrag_back(2,i,iset) if(me.eq.king.or..not.out1file) & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset), & wfrag_back(2,i,iset),qin_back(3,i,iset), & wfrag_back(3,i,iset),qin_back(3,i,iset), & ifrag_back(1,i,iset),ifrag_back(2,i,iset) enddo else do i=1,nfrag_back read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset), & wfrag_back(3,i,iset), & ifrag_back(1,i,iset),ifrag_back(2,i,iset) if(me.eq.king.or..not.out1file) & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset), & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset) enddo endif enddo return end C--------------------------------------------------------------------------- subroutine read_afminp implicit none 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 integer i c 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) c print *,'initdist',distafminit return end c------------------------------------------------------------------------------- subroutine read_saxs_constr implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.SAXS' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' double precision cm(3),cnorm integer i,j c read(inp,*) nsaxs write (iout,*) "Calling read_saxs nsaxs",nsaxs call flush(iout) if (saxs_mode.eq.0) then c SAXS distance distribution do i=1,nsaxs read(inp,*) distsaxs(i),Psaxs(i) enddo Cnorm = 0.0d0 do i=1,nsaxs Cnorm = Cnorm + Psaxs(i) enddo write (iout,*) "Cnorm",Cnorm do i=1,nsaxs Psaxs(i)=Psaxs(i)/Cnorm enddo write (iout,*) "Normalized distance distribution from SAXS" do i=1,nsaxs write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i) enddo Wsaxs0=0.0d0 do i=1,nsaxs Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i)) enddo write (iout,*) "Wsaxs0",Wsaxs0 else c SAXS "spheres". do i=1,nsaxs read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3) enddo do j=1,3 cm(j)=0.0d0 enddo do i=1,nsaxs do j=1,3 cm(j)=cm(j)+Csaxs(j,i) enddo enddo do j=1,3 cm(j)=cm(j)/nsaxs enddo do i=1,nsaxs do j=1,3 Csaxs(j,i)=Csaxs(j,i)-cm(j) enddo enddo write (iout,*) "SAXS sphere coordinates" do i=1,nsaxs write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3) enddo endif return end c------------------------------------------------------------------------------- subroutine read_dist_constr implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' include 'COMMON.INTERACT' integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000) double precision wfrag_(100),wpair_(1000) double precision ddjk,dist,dist_cut,fordepthmax character*5000 controlcard logical normalize,next integer restr_type double precision scal_bfac double precision xlink(4,0:4) / c a b c sigma & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS c print *, "WCHODZE" write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) restr_on_coord=.false. next=.true. npeak=0 ipeak=0 nhpb_peak=0 DO WHILE (next) call card_concat(controlcard) next = index(controlcard,"NEXT").gt.0 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist) write (iout,*) "restr_type",restr_type call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) call readi(controlcard,"NDIST",ndist_,0) call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0) if (restr_type.eq.10) & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0) if (restr_type.eq.12) & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0) call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) normalize = index(controlcard,"NORMALIZE").gt.0 write (iout,*) "WBOLTZD",wboltzd write (iout,*) "SCAL_PEAK",scal_peak write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ write (iout,*) "IFRAG" do i=1,nfrag_ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) enddo write (iout,*) "IPAIR" do i=1,npair_ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) enddo if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) & write (iout,*) & "Distance restraints as generated from reference structure" do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup if (ifrag_(2,i).gt.nstart_sup+nsup-1) & ifrag_(2,i)=nstart_sup+nsup-1 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) c call flush(iout) if (wfrag_(i).eq.0.0d0) cycle do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo enddo do i=1,npair_ if (wpair_(i).eq.0.0d0) cycle ii = ipair_(1,i) jj = ipair_(2,i) if (ii.gt.jj) then itemp=ii ii=jj jj=itemp endif do j=ifrag_(1,ii),ifrag_(2,ii) do k=ifrag_(1,jj),ifrag_(2,jj) ddjk=dist(j,k) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wpair_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wpair_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo enddo c print *,ndist_ write (iout,*) "Distance restraints as read from input" do i=1,ndist_ if (restr_type.eq.12) then read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), & fordepth_peak(nhpb_peak+1),npeak c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), c & fordepth_peak(nhpb_peak+1),npeak if (forcon_peak(nhpb_peak+1).le.0.0d0.or. & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle nhpb_peak=nhpb_peak+1 irestr_type_peak(nhpb_peak)=12 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i ipeak(2,npeak)=i #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) #else write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) #endif if (ibecarb_peak(nhpb_peak).eq.3) then jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres else if (ibecarb_peak(nhpb_peak).eq.2) then ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres else if (ibecarb_peak(nhpb_peak).eq.1) then ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres endif else if (restr_type.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle nhpb=nhpb+1 irestr_type(nhpb)=11 #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) #else write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) #endif c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif else if (restr_type.eq.10) then c Cross-lonk Markov-like potential call card_concat(controlcard) call readi(controlcard,"ILINK",ihpb(nhpb+1),0) call readi(controlcard,"JLINK",jhpb(nhpb+1),0) ibecarb(nhpb+1)=0 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle if (index(controlcard,"ZL").gt.0) then link_type=1 else if (index(controlcard,"ADH").gt.0) then link_type=2 else if (index(controlcard,"PDH").gt.0) then link_type=3 else if (index(controlcard,"DSS").gt.0) then link_type=4 else link_type=0 endif call reada(controlcard,"AXLINK",dhpb(nhpb+1), & xlink(1,link_type)) call reada(controlcard,"BXLINK",dhpb1(nhpb+1), & xlink(2,link_type)) call reada(controlcard,"CXLINK",fordepth(nhpb+1), & xlink(3,link_type)) call reada(controlcard,"SIGMA",forcon(nhpb+1), & xlink(4,link_type)) call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0) c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1), c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) if (forcon(nhpb+1).le.0.0d0 .or. & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle nhpb=nhpb+1 irestr_type(nhpb)=10 if (ibecarb(nhpb).eq.3) then jhpb(nhpb)=jhpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) #else write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) #endif else C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 if (dhpb1(nhpb).eq.0.0d0) then irestr_type(nhpb)=1 else irestr_type(nhpb)=2 endif if (ibecarb(nhpb).eq.3) then jhpb(nhpb)=jhpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif if (dhpb(nhpb).eq.0.0d0) & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) endif #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #endif endif C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) C if (forcon(nhpb+1).gt.0.0d0) then C nhpb=nhpb+1 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) enddo if (restr_type.eq.4) then write (iout,*) "The BFAC array" do i=nnt,nct write (iout,'(i5,f10.5)') i,bfac(i) enddo do i=nnt,nct if (itype(i).eq.ntyp1) cycle do j=nnt,i-1 if (itype(j).eq.ntyp1) cycle if (itype(i).eq.10) then iiend=0 else iiend=1 endif if (itype(j).eq.10) then jjend=0 else jjend=1 endif kk=0 do ii=0,iiend do jj=0,jjend nhpb=nhpb+1 irestr_type(nhpb)=1 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2) irestr_type(nhpb)=1 ibecarb(nhpb)=kk if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb) ihpb(nhpb)=i+nres*ii jhpb(nhpb)=j+nres*jj dhpb(nhpb)=dist(i+nres*ii,j+nres*jj) #ifdef MPI if (.not.out1file .or. me.eq.king) then write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) endif #else write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) #endif kk=kk+1 enddo enddo enddo enddo endif if (restr_type.eq.5) then restr_on_coord=.true. do i=nnt,nct if (itype(i).eq.ntyp1) cycle bfac(i)=(scal_bfac/bfac(i))**2 enddo endif ENDDO ! next fordepthmax=0.0d0 if (normalize) then do i=nss+1,nhpb if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) & fordepthmax=fordepth(i) enddo do i=nss+1,nhpb if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax enddo endif return end c------------------------------------------------------------------------------- subroutine read_constr_homology implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.HOMOLOGY' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.MD' include 'COMMON.QRESTR' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.NAMES' c c For new homol impl c include 'COMMON.VAR' c c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, c & dist_cut c common /przechowalnia/ odl_temp(maxres,maxres,max_template), c & sigma_odl_temp(maxres,maxres,max_template) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec, & ik,iistart integer ilen external ilen logical liiflag,lfirst integer i01,i10 c c FP - Nov. 2014 Temporary specifications for new vars c double precision rescore_tmp,x12,y12,z12,rescore2_tmp, & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 double precision, dimension (max_template,maxres) :: rescore3 double precision distal character*24 pdbfile,tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints c not using pre-computed ones stored in files model_ki_{dist,angle} c FP (Nov., 2014) c ----------------------------------------------------------------- c c c Alternative: reading from input call card_concat(controlcard) call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0) call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0) call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) if(.not.read2sigma.and.start_from_model) then if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA' start_from_model=.false. endif if(start_from_model .and. (me.eq.king .or. .not. out1file)) & write(iout,*) 'START_FROM_MODELS is ON' if(start_from_model .and. rest) then if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) 'START_FROM_MODELS is OFF' write(iout,*) 'remove restart keyword from input' endif endif if (homol_nset.gt.1)then call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "iset homology_weight " do i=1,homol_nset write(iout,*) i,waga_homology(i) enddo endif iset=mod(kolor,homol_nset)+1 else iset=1 waga_homology(1)=1.0 endif cd write (iout,*) "nnt",nnt," nct",nct cd call flush(iout) lim_odl=0 lim_dih=0 c c write(iout,*) 'nnt=',nnt,'nct=',nct c do i = nnt,nct do k=1,constr_homology idomain(k,i)=0 enddo enddo ii=0 do i = nnt,nct-2 do j=i+2,nct ii=ii+1 ii_in_use(ii)=0 enddo enddo if (read_homol_frag) then call read_klapaucjusz else do k=1,constr_homology read(inp,'(a)') pdbfile if(me.eq.king .or. .not. out1file) & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', & pdbfile(:ilen(pdbfile)) stop 34 continue c print *,'Begin reading pdb data' c c Files containing res sim or local scores (former containing sigmas) c write(kic2,'(bz,i2.2)') k tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. if (read2sigma) then call readpdb_template(k) else call readpdb endif c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) do i=1,2*nres do j=1,3 c(j,i)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo c c From read_dist_constr (commented out 25/11/2014 <-> res sim) c c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore open (ientin,file=tpl_k_rescore,status='old') if (nnt.gt.1) rescore(k,1)=0.0d0 do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp rescore3(k,i_tmp)=rescore3_tmp if (.not. out1file .or. me.eq.king) & write(iout,'(a7,i5,3f10.5,i5)') "rescore", & i_tmp,rescore2_tmp,rescore_tmp, & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif enddo 1401 continue close (ientin) if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", c & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif c write (iout,*) "Distance restraints set" c call flush(iout) c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_dih,status='old') c do irec=1,maxres-3 ! loop for reading sigma_dih c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_dih(k,i+nnt-1) c enddo c1402 continue c close (ientin) do i = nnt+3,nct if (idomain(k,i).eq.0) then sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) ! right? c read (ientin,*) sigma_dih(k,i) ! original variant c write (iout,*) "dih(",k,i,") =",dih(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2), c & "rescore(",k,i-3,") =",rescore(k,i-3) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 c write (iout,*) "Raw sigmas for dihedral angle restraints" c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2)*rescore(k,i-3) ! right expression ? c Instead of res sim other local measure of b/b str reliability possible if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 endif c write (iout,*) "Dihedral angle restraints set" c call flush(iout) if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_theta(k,i+nnt-1) c enddo c1403 continue c close (ientin) do i = nnt+2,nct ! right? without parallel. c do i = i=1,nres ! alternative for bounds acc to readpdb? c do i=ithet_start,ithet_end ! with FG parallel. if (idomain(k,i).eq.0) then sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2) c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif c write (iout,*) "Angle restraints set" c call flush(iout) if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_d(k,i+nnt-1) c enddo c1404 continue do i = nnt,nct ! right? without parallel. c do i=2,nres-1 ! alternative for bounds acc to readpdb? c do i=loc_start,loc_end ! with FG parallel. if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i) sigma_d(k,i)=rescore3(k,i) ! right expression ? if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) c read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo c write (iout,*) "SC restraints set" c call flush(iout) c c remove distance restraints not used in any model from the list c shift data in all arrays c c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. lfirst=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 c & .and. distal.le.dist2_cut ) then c write (iout,*) "i",i," j",j," ii",ii c call flush(iout) if (ii_in_use(ii).eq.0.and.liiflag.or. & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then liiflag=.false. i10=ii if (lfirst) then lfirst=.false. iistart=ii else if(i10.eq.lim_odl) i10=i10+1 do ki=0,i10-i01-1 ires_homo(iistart+ki)=ires_homo(ki+i01) jres_homo(iistart+ki)=jres_homo(ki+i01) ii_in_use(iistart+ki)=ii_in_use(ki+i01) do k=1,constr_homology odl(k,iistart+ki)=odl(k,ki+i01) sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) l_homo(k,iistart+ki)=l_homo(k,ki+i01) enddo enddo iistart=iistart+i10-i01 endif endif if (ii_in_use(ii).ne.0.and..not.liiflag) then i01=ii liiflag=.true. endif enddo enddo lim_odl=iistart-1 endif c write (iout,*) "Removing distances completed" c call flush(iout) endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition c write (iout,*) "After homology_partition" c call flush(iout) if (constr_homology.gt.0) call init_int_table c write (iout,*) "After init_int_table" c call flush(iout) c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end c c Print restraints c if (.not.out_template_restr) return cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" do ii=1,lim_odl write(iout,'(3i5,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii), & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" do i=nnt+3,nct write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" do i=nnt+2,nct write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) enddo endif c ----------------------------------------------------------------- return end c---------------------------------------------------------------------- #ifdef WINIFL subroutine flush(iu) return end #endif #ifdef AIX subroutine flush(iu) call flush_(iu) return end #endif c------------------------------------------------------------------------------ subroutine copy_to_tmp(source) implicit none include "DIMENSIONS" include "COMMON.IOUNITS" character*(*) source character* 256 tmpfile integer ilen external ilen logical ex tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source)) inquire(file=tmpfile,exist=ex) if (ex) then write (*,*) "Copying ",tmpfile(:ilen(tmpfile)), & " to temporary directory..." write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir) endif return end c------------------------------------------------------------------------------ subroutine move_from_tmp(source) implicit none include "DIMENSIONS" include "COMMON.IOUNITS" character*(*) source integer ilen external ilen write (*,*) "Moving ",source(:ilen(source)), & " from temporary directory to working directory" write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir call system("/bin/mv "//source(:ilen(source))//" "//curdir) return end c------------------------------------------------------------------------------ subroutine random_init(seed) C C Initialize random number generator C implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' logical OKRandom, prng_restart real*8 r1 integer iseed_array(4) integer error_msg,ierr #endif include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.THREAD' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.MCM' include 'COMMON.MAP' include 'COMMON.HEADER' include 'COMMON.CSA' include 'COMMON.CHAIN' include 'COMMON.MUCA' include 'COMMON.MD' include 'COMMON.FFIELD' include 'COMMON.SETUP' integer i,iseed double precision seed,ran_number iseed=-dint(dabs(seed)) if (iseed.eq.0) then write (iout,'(/80(1h*)/20x,a/80(1h*))') & 'Random seed undefined. The program will stop.' write (*,'(/80(1h*)/20x,a/80(1h*))') & 'Random seed undefined. The program will stop.' #ifdef MPI call mpi_finalize(mpi_comm_world,ierr) #endif stop 'Bad random seed.' endif #ifdef MPI if (fg_rank.eq.0) then seed=seed*(me+1)+1 #ifdef AMD64 if(me.eq.king) & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed OKRandom = prng_restart(me,iseed) #else do i=1,4 tmp=65536.0d0**(4-i) iseed_array(i) = dint(seed/tmp) seed=seed-iseed_array(i)*tmp enddo if(me.eq.king) & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ', & (iseed_array(i),i=1,4) write (*,*) 'MPI: node= ',me, ' iseed(4)= ', & (iseed_array(i),i=1,4) OKRandom = prng_restart(me,iseed_array) #endif if (OKRandom) then c r1 = prng_next(me) r1=ran_number(0.0D0,1.0D0) if(me.eq.king) & write (iout,*) 'ran_num',r1 if (r1.lt.0.0d0) OKRandom=.false. endif if (.not.OKRandom) then write (iout,*) 'PRNG IS NOT WORKING!!!' print *,'PRNG IS NOT WORKING!!!' if (me.eq.0) then call flush(iout) call mpi_abort(mpi_comm_world,error_msg,ierr) stop else write (iout,*) 'too many processors for parallel prng' write (*,*) 'too many processors for parallel prng' call flush(iout) stop endif endif endif #else call vrndst(iseed) write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0) #endif return end c---------------------------------------------------------------------- subroutine read_klapaucjusz implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.HOMOLOGY' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.MD' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.NAMES' character*256 fragfile integer ninclust(maxclust),inclust(max_template,maxclust), & nresclust(maxclust),iresclust(maxres,maxclust),nclust character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp, & ik,ll,ii,kk,iistart,iishift,lim_xx double precision distal logical lprn /.true./ integer ilen external ilen logical liiflag c c double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 character*24 pdbfile,tpl_k_rescore c c For new homol impl c include 'COMMON.VAR' c call getenv("FRAGFILE",fragfile) open(ientin,file=fragfile,status="old",err=10) read(ientin,*) constr_homology,nclust l_homo = .false. sigma_theta=0.0 sigma_d=0.0 sigma_dih=0.0 c Read pdb files do k=1,constr_homology read(ientin,'(a)') pdbfile write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', & pdbfile(:ilen(pdbfile)) stop 34 continue unres_pdb=.false. call readpdb_template(k) do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 enddo enddo c Read clusters do i=1,nclust read(ientin,*) ninclust(i),nresclust(i) read(ientin,*) (inclust(k,i),k=1,ninclust(i)) read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) enddo c c Loop over clusters c do l=1,nclust do ll = 1,ninclust(l) k = inclust(ll,l) do i=1,nres idomain(k,i)=0 enddo do i=1,nresclust(l) if (nnt.gt.1) then idomain(k,iresclust(i,l)+1) = 1 else idomain(k,iresclust(i,l)) = 1 endif enddo c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", c & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 c l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then do i = nnt+3,nct if (idomain(k,i).eq.0) then c sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), c & " sigma_dihed",sigma_dih(k,i) if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) enddo lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then do i = nnt+2,nct if (idomain(k,i).eq.0) then c sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) enddo endif if (waga_d.gt.0.0d0) then do i = nnt,nct if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then c sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) sigma_d(k,i)=rescore(k,i) if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 enddo endif enddo ! l enddo ! ll c c remove distance restraints not used in any model from the list c shift data in all arrays c if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 if (ii_in_use(ii).eq.0.and.liiflag) then liiflag=.false. iistart=ii endif if (ii_in_use(ii).ne.0.and..not.liiflag.or. & .not.liiflag.and.ii.eq.lim_odl) then if (ii.eq.lim_odl) then iishift=ii-iistart+1 else iishift=ii-iistart endif liiflag=.true. do ki=iistart,lim_odl-iishift ires_homo(ki)=ires_homo(ki+iishift) jres_homo(ki)=jres_homo(ki+iishift) ii_in_use(ki)=ii_in_use(ki+iishift) do k=1,constr_homology odl(k,ki)=odl(k,ki+iishift) sigma_odl(k,ki)=sigma_odl(k,ki+iishift) l_homo(k,ki)=l_homo(k,ki+iishift) enddo enddo ii=ii-iishift lim_odl=lim_odl-iishift endif enddo enddo endif return 10 stop "Error in fragment file" end