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) 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 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' 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,0) 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) C Juyong call reada(weightcard,'WDFAD',wdfa_dist,0.0d0) call reada(weightcard,'WDFAT',wdfa_tor,0.0d0) call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) C 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 C JUYONG weights(24)=wdfa_dist weights(25)=wdfa_tor weights(26)=wdfa_nei weights(27)=wdfa_beta C 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, & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 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)'/ & 'WDFA_D= ',f10.6,' (DFA, distance)' / & 'WDFA_T= ',f10.6,' (DFA, torsional)' / & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / & 'WDFA_B= ',f10.6,' (DFA, beta formation)') 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 c 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, & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 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)'/ & 'WDFA_D= ',f10.6,' (DFA, distance)' / & 'WDFA_T= ',f10.6,' (DFA, torsional)' / & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / & 'WDFA_B= ',f10.6,' (DFA, beta formation)') 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 call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then C Following 2 lines for diagnostics; comment out if not needed write (iout,*) "Before sideadd" call intout 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 C 10/03/12 Adam: Recalculate coordinates with new side chain positions call chainbuild C Following 2 lines for diagnostics; comment out if not needed write (iout,*) "After sideadd" call intout 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 C Juyong:READ init_vars C Initialize variables! C Juyong:READ read_info C READ fragment information!! C both routines should be in dfa.F file!! 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 C C 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 call contact(.true.,ncont_ref,icont_ref,co) endif 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 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_dist.gt.0) then call read_dist_constr endif if (nhpb.gt.0) call hpb_partition c write (iout,*) "After read_dist_constr nhpb",nhpb c call flush(iout) 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. 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 c---------------------------------------------------------------------------- subroutine csaread implicit real*8 (a-h,o-z) 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,-300000.0d0) call readi(mcmcard,'ICMAX',icmax,1) call readi(mcmcard,'NBANKM',nbankm,400) call readi(mcmcard,'IUCUT',iucut,2) 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 tm_score=(index(mcmcard,'TMSCORE').gt.0) if (tm_score) write (iout,*) "Using TM_Score instead of DIFF", & " for torsional angles" 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) 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 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 if (.not.refstr .and. nfrag.gt.0) then write (iout,*) & "ERROR: no reference structure to compute distance restraints" write (iout,*) & "Restraints must be specified explicitly (NDIST=number)" stop endif if (nfrag.lt.2 .and. npair.gt.0) then write (iout,*) "ERROR: Less than 2 fragments specified", & " but distance restraints between pairs requested" stop endif 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) c 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),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 if (ibecarb(i).gt.0) then ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres endif if (dhpb(nhpb).eq.0.0d0) & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) endif enddo #ifdef MPI if (.not.out1file .or. me.eq.king) then #endif do i=1,nhpb write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) enddo call flush(iout) #ifdef MPI endif #endif 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 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) write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0) #endif return end