subroutine readrtns implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' logical file_exist C Read force-field parameters except weights call parmread C Read job setup parameters call read_control C Read control parameters for energy minimzation if required if (minim) call read_minim C Read MCM control parameters if required c if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread C Read MD control parameters if reqjuired c if (modecalc.eq.12) call read_MDpar C Read MREMD control parameters if required c if (modecalc.eq.14) then c call read_MDpar c call read_REMDpar c endif C Read MUCA control parameters if required c if (lmuca) call read_muca C Read CSA control parameters if required (from fort.40 if exists C otherwise from general input file) c if (modecalc.eq.8) then c inquire (file="fort.40",exist=file_exist) c if (.not.file_exist) call csaread c endif cfmc if (modecalc.eq.10) call mcmfread 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) &write (iout,'(a,i5,a)') "The following",nhpb-nss, & " distance constraints have been imposed" do i=nss+1,nhpb write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i) enddo #ifdef MPI endif #endif c print *,"Processor",myrank," leaves READRTNS" return end C------------------------------------------------------------------------------- subroutine read_control C C Read contorl data C implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MP include 'mpif.h' logical OKRandom, prng_restart real*8 r1 #endif include 'COMMON.IOUNITS' include 'COMMON.TIME1' c include 'COMMON.THREAD' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.MCM' c include 'COMMON.MAP' include 'COMMON.HEADER' c include 'COMMON.CSA' include 'COMMON.CHAIN' c include 'COMMON.MUCA' c include 'COMMON.MD' include 'COMMON.FFIELD' include 'COMMON.SETUP' COMMON /MACHSW/ KDIAG,ICORFL,IXDR character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ character*80 ucase character*320 controlcard 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 call readi(controlcard,'CONSTR_DIST',constr_dist,0) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes call reada(controlcard,'RMSDBC',rmsdbc,3.0D0) call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0) call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0) call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0) call reada(controlcard,'DRMS',drms,0.1D0) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max write (iout,'(a,f10.1)')'DRMS = ',drms write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm write (iout,'(a,f10.1)') 'Time limit (min):',timlim endif call readi(controlcard,'NZ_START',nz_start,0) call readi(controlcard,'NZ_END',nz_end,0) call readi(controlcard,'IZ_SC',iz_sc,0) timlim=60.0D0*timlim safety = 60.0d0*safety timem=timlim modecalc=0 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,'NOSEARCHSC').gt.0) searchsc=.not.searchsc sideadd=(index(controlcard,'SIDEADD').gt.0) energy_dec=(index(controlcard,'ENERGY_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) call readi(controlcard,'IPRINT',iprint,0) 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,1) 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 (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 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 if(me.eq.king.or..not.out1file) & write (iout,'(2a)') diagmeth(kdiag), & ' routine used to diagonalize matrices.' return end c------------------------------------------------------------------------------ subroutine molread C C Read molecular data. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer error_msg #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' c include 'COMMON.DBASE' c include 'COMMON.THREAD' include 'COMMON.CONTACTS' include 'COMMON.TORCNSTR' include 'COMMON.TIME1' include 'COMMON.BOUNDS' c include 'COMMON.MD' c include 'COMMON.REMD' include 'COMMON.SETUP' character*4 sequence(maxres) integer rescode double precision x(maxvar) character*256 pdbfile character*320 weightcard character*80 weightcard_t,ucase dimension itype_pdb(maxres) common /pizda/ itype_pdb logical seq_comp,fail double precision energia(0:n_ene) integer ilen external ilen C C Body C C Read weights of the subsequent energy terms. call card_concat(weightcard) call reada(weightcard,'WLONG',wlong,1.0D0) call reada(weightcard,'WSC',wsc,wlong) call reada(weightcard,'WSCP',wscp,wlong) call reada(weightcard,'WELEC',welec,1.0D0) call reada(weightcard,'WVDWPP',wvdwpp,welec) call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) call reada(weightcard,'WCORR4',wcorr4,0.0D0) call reada(weightcard,'WCORR5',wcorr5,0.0D0) call reada(weightcard,'WCORR6',wcorr6,0.0D0) call reada(weightcard,'WTURN3',wturn3,1.0D0) call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp weights(3)=welec weights(4)=wcorr weights(5)=wcorr5 weights(6)=wcorr6 weights(7)=wel_loc weights(8)=wturn3 weights(9)=wturn4 weights(10)=wturn6 weights(11)=wang weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d weights(15)=wstrain weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 weights(21)=wsccor if(me.eq.king.or..not.out1file) & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, & wturn4,wturn6 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)') if(me.eq.king.or..not.out1file)then if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' write (iout,'(2(a,f5.3/))') & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr, & 'Range of quenching the correlation terms:',2*delt_corr else if (wcorr.gt.0.0d0) then write (iout,'(/2a/)') 'Hydrogen-bonding correlation ', & 'between contact pairs of peptide groups' endif write (iout,'(a,f8.3)') & 'Scaling factor of 1,4 SC-p interactions:',scal14 write (iout,'(a,f8.3)') & 'General scaling factor of SC-p interactions:',scalscp endif r0_corr=cutoff_corr-delt_corr do i=1,20 aad(i,1)=scalscp*aad(i,1) aad(i,2)=scalscp*aad(i,2) bad(i,1)=scalscp*bad(i,1) bad(i,2)=scalscp*bad(i,2) enddo call rescale_weights(t_bath) if(me.eq.king.or..not.out1file) & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, & wturn4,wturn6 22 format (/'Energy-term weights (scaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)') if(me.eq.king.or..not.out1file) & write (iout,*) "Reference temperature for weights calculation:", & temp0 call reada(weightcard,"D0CM",d0cm,3.78d0) call reada(weightcard,"AKCM",akcm,15.1d0) call reada(weightcard,"AKTH",akth,11.0d0) call reada(weightcard,"AKCT",akct,12.0d0) call reada(weightcard,"V1SS",v1ss,-1.08d0) call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) if(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, & " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss write (iout,*) "EBR",ebr print *,'indpdb=',indpdb,' pdbref=',pdbref endif 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 c print *,'Begin reading pdb data' call readpdb c print *,'Finished reading pdb data' 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 c 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) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) c 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 C Assign initial virtual bond lengths do i=2,nres vbld(i)=vbl vbld_inv(i)=vblinv enddo do i=2,nres-1 vbld(i+nres)=dsc(itype(i)) vbld_inv(i+nres)=dsc_inv(itype(i)) c write (iout,*) "i",i," itype",itype(i), c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres) enddo endif c print *,nres c print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR if (itype(i).eq.21 .or. itype(i+1).eq.21) then #else if (itype(i).eq.21) then #endif itel(i)=0 #ifdef PROCOR else if (itype(i+1).ne.20) then #else else if (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 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 if (ndih_constr.gt.0) then read (inp,*) ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) enddo endif do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo if(me.eq.king.or..not.out1file) & 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 endif nnt=1 #ifdef MPI if (me.eq.king) then #endif write (iout,'(a)') 'Boundaries in phi 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 nct=nres cd print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) nct=nct-1 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 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) 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 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 c call contact(.true.,ncont_ref,icont_ref,co) endif c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr c write (iout,*) "After read_dist_constr nhpb",nhpb call hpb_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 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) 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_norm(j,i)*vbld_inv(i+1) enddo enddo do i=nnt,nct if (itype(i).ne.10) then do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres) enddo endif enddo return else call read_angles(inp,*36) 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 enddo else if(me.eq.king.or..not.out1file) & write (iout,'(a)') 'Random-generated initial geometry.' #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 c 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 do itrial=1,100 itmp=1 c call gen_rand_conf(itmp,*31) goto 40 31 write (iout,*) 'Failed to generate random conformation', & ', itrial=',itrial write (*,*) 'Failed to generate random conformation', & ', itrial=',itrial enddo write (iout,'(a,i3,a)') 'Processor:',me, & ' error in generating random conformation.' write (*,'(a,i3,a)') 'Processor:',me, & ' 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) 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) if (me.eq.king.or..not.out1file) & 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 c 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 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 real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #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' c include 'COMMON.DBASE' c include 'COMMON.THREAD' include 'COMMON.TIME1' include 'COMMON.SETUP' 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(iss(i)),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ',restyp(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) 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 dhpb(i)=dbr 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 real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.LOCAL' include 'COMMON.INTERACT' 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) 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 setup_var implicit real*8 (a-h,o-z) 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' c include 'COMMON.DBASE' c include 'COMMON.THREAD' include 'COMMON.TIME1' 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) 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 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1) return end c---------------------------------------------------------------------------- subroutine gen_dist_constr C Generate CA distance constraints. implicit real*8 (a-h,o-z) 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' c include 'COMMON.DBASE' c include 'COMMON.THREAD' include 'COMMON.TIME1' dimension itype_pdb(maxres) common /pizda/ itype_pdb 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 subroutine read_minim implicit real*8 (a-h,o-z) 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) 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 real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' 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 real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' character*16 form,nodename integer nodelen #endif include 'COMMON.SETUP' include 'COMMON.IOUNITS' c include 'COMMON.MD' include 'COMMON.CONTROL' integer lenpre,lenpot,ilen,lentmp external ilen character*3 out1file_text,ucase character*3 ll external ucase 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) #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared) #endif call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',readonly,shared) #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared) #endif 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) #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') c print *,"Processor",myrank," opened file ITHEP" #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) open (ithep_pdb,file=thetname_pdb,status='old',action='read') #endif call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',action='read') c print *,"Processor",myrank," opened file IROTAM" #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif 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') 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') #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) #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) print *,"thetname_pdb ",thetname_pdb open (ithep_pdb,file=thetname_pdb,status='old',readonly) print *,ithep_pdb," opened" #endif call getenv_loc('ROTPAR',rotname) open (irotam,file=rotname,status='old',readonly) #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',readonly) #endif 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) 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) #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' c if(usampl) then c qname=prefix(:lenpre)//'_'//pot(:lenpot)// c & liczba(:ilen(liczba))//'.const' c 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' c if(usampl) then c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const' c endif #endif #if defined(AIX) || defined(PGI) 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',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 parameter file : ", & thetname(:ilen(thetname)) write (iout,*) "Rotamer parameter file : ", & rotname(:ilen(rotname)) 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 subroutine read_dist_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard c write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) call card_concat(controlcard) 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 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) c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ c write (iout,*) "IFRAG" c do i=1,nfrag_ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) c enddo c write (iout,*) "IPAIR" c do i=1,npair_ c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) c enddo call flush(iout) 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) call flush(iout) if (wfrag_(i).gt.0.0d0) then do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) write (iout,*) "j",j," k",k ddjk=dist(j,k) if (constr_dist.eq.1) then nhpb=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 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) endif else nhpb=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.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo endif enddo do i=1,npair_ if (wpair_(i).gt.0.0d0) then 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) nhpb=nhpb+1 ihpb(nhpb)=j jhpb(nhpb)=k forcon(nhpb)=wpair_(i) dhpb(nhpb)=dist(j,k) #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo endif enddo do i=1,ndist_ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif endif enddo call flush(iout) 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) 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) 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 real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef AMD64 integer*8 iseedi8 #endif #ifdef MPI include 'mpif.h' logical OKRandom, prng_restart real*8 r1 integer iseed_array(4) #endif include 'COMMON.IOUNITS' include 'COMMON.TIME1' c include 'COMMON.THREAD' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.MCM' c include 'COMMON.MAP' include 'COMMON.HEADER' c include 'COMMON.CSA' include 'COMMON.CHAIN' c include 'COMMON.MUCA' c include 'COMMON.MD' include 'COMMON.FFIELD' include 'COMMON.SETUP' 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 iseedi8=dint(seed) if(me.eq.king .or. .not. out1file) & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8 OKRandom = prng_restart(me,iseedi8) #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 .or. .not. out1file) & 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=ran_number(0.0D0,1.0D0) if(me.eq.king .or. .not. out1file) & 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) c write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0) #endif return end