2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
16 C Read control parameters for energy minimzation if required
17 if (minim) call read_minim
18 C Read molecule information, molecule geometry, energy-term weights, and
19 C restraints if requested
21 C Print restraint information
23 if (.not. out1file .or. me.eq.king) then
26 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
27 & " distance constraints have been imposed"
29 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
30 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
35 c print *,"Processor",myrank," leaves READRTNS"
38 C-------------------------------------------------------------------------------
39 subroutine read_control
43 implicit real*8 (a-h,o-z)
47 logical OKRandom, prng_restart
50 include 'COMMON.IOUNITS'
51 include 'COMMON.TIME1'
52 include 'COMMON.THREAD'
53 include 'COMMON.SBRIDGE'
54 include 'COMMON.CONTROL'
57 include 'COMMON.HEADER'
58 csa include 'COMMON.CSA'
59 include 'COMMON.CHAIN'
62 include 'COMMON.FFIELD'
63 include 'COMMON.SETUP'
64 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
65 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
67 character*320 controlcard
72 read (INP,'(a)') titel
73 call card_concat(controlcard)
74 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
75 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
76 call reada(controlcard,'SEED',seed,0.0D0)
77 call random_init(seed)
78 C Set up the time limit (caution! The time must be input in minutes!)
79 read_cart=index(controlcard,'READ_CART').gt.0
80 catrace=index(controlcard,'CATRACE').gt.0
81 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
82 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
83 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
84 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
85 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
86 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
87 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
88 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
89 call reada(controlcard,'DRMS',drms,0.1D0)
90 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
91 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
92 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
93 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
94 write (iout,'(a,f10.1)')'DRMS = ',drms
95 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
96 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
98 call readi(controlcard,'NZ_START',nz_start,0)
99 call readi(controlcard,'NZ_END',nz_end,0)
100 call readi(controlcard,'IZ_SC',iz_sc,0)
102 safety = 60.0d0*safety
105 call reada(controlcard,"T_BATH",t_bath,300.0d0)
106 minim=(index(controlcard,'MINIMIZE').gt.0)
107 dccart=(index(controlcard,'CART').gt.0)
108 overlapsc=(index(controlcard,'OVERLAP').gt.0)
109 overlapsc=.not.overlapsc
110 searchsc=(index(controlcard,'SEARCHSC').gt.0)
111 sideadd=(index(controlcard,'SIDEADD').gt.0) .or. catrace
112 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
113 outpdb=(index(controlcard,'PDBOUT').gt.0)
114 outmol2=(index(controlcard,'MOL2OUT').gt.0)
115 pdbref=(index(controlcard,'PDBREF').gt.0)
116 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
117 indpdb=index(controlcard,'PDBSTART')
118 extconf=(index(controlcard,'EXTCONF').gt.0)
119 call readi(controlcard,'IPRINT',iprint,0)
120 call readi(controlcard,'MAXGEN',maxgen,10000)
121 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
122 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
123 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
124 & write (iout,*) "RESCALE_MODE",rescale_mode
125 if (index(controlcard,'REGULAR').gt.0.0D0) then
126 call reada(controlcard,'WEIDIS',weidis,0.1D0)
131 if (index(controlcard,'CHECKGRAD').gt.0) then
133 if (index(controlcard,'CART').gt.0) then
135 elseif (index(controlcard,'CARINT').gt.0) then
140 elseif (index(controlcard,'THREAD').gt.0) then
142 call readi(controlcard,'THREAD',nthread,0)
143 if (nthread.gt.0) then
144 call reada(controlcard,'WEIDIS',weidis,0.1D0)
147 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
148 stop 'Error termination in Read_Control.'
150 else if (index(controlcard,'MCMA').gt.0) then
152 else if (index(controlcard,'MCEE').gt.0) then
154 else if (index(controlcard,'MULTCONF').gt.0) then
156 else if (index(controlcard,'MAP').gt.0) then
158 call readi(controlcard,'MAP',nmap,0)
159 else if (index(controlcard,'CSA').gt.0) then
160 write(*,*) "CSA not supported in this version"
162 else if (index(controlcard,'SOFTREG').gt.0) then
164 else if (index(controlcard,'CHECK_BOND').gt.0) then
166 else if (index(controlcard,'TEST').gt.0) then
170 lmuca=index(controlcard,'MUCA').gt.0
171 call readi(controlcard,'MUCADYN',mucadyn,0)
172 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
173 if (lmuca .and. (me.eq.king .or. .not.out1file ))
175 write (iout,*) 'MUCADYN=',mucadyn
176 write (iout,*) 'MUCASMOOTH=',muca_smooth
179 iscode=index(controlcard,'ONE_LETTER')
180 indphi=index(controlcard,'PHI')
181 indback=index(controlcard,'BACK')
182 iranconf=index(controlcard,'RAND_CONF')
183 i2ndstr=index(controlcard,'USE_SEC_PRED')
184 gradout=index(controlcard,'GRADOUT').gt.0
185 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
187 if(me.eq.king.or..not.out1file)
188 & write (iout,'(2a)') diagmeth(kdiag),
189 & ' routine used to diagonalize matrices.'
192 c--------------------------------------------------------------------------
195 C Read molecular data.
197 implicit real*8 (a-h,o-z)
203 include 'COMMON.IOUNITS'
206 include 'COMMON.INTERACT'
207 include 'COMMON.LOCAL'
208 include 'COMMON.NAMES'
209 include 'COMMON.CHAIN'
210 include 'COMMON.FFIELD'
211 include 'COMMON.SBRIDGE'
212 include 'COMMON.HEADER'
213 include 'COMMON.CONTROL'
214 include 'COMMON.DBASE'
215 include 'COMMON.THREAD'
216 include 'COMMON.CONTACTS'
217 include 'COMMON.TORCNSTR'
218 include 'COMMON.TIME1'
219 include 'COMMON.BOUNDS'
221 include 'COMMON.REMD'
222 include 'COMMON.SETUP'
223 character*4 sequence(maxres)
225 double precision x(maxvar)
226 character*320 weightcard
227 character*80 weightcard_t,ucase
228 dimension itype_pdb(maxres)
229 common /pizda/ itype_pdb
230 logical seq_comp,fail
231 double precision energia(0:n_ene)
237 C Read weights of the subsequent energy terms.
250 if(me.eq.king.or..not.out1file) then
251 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
252 write (iout,*) 'Current weights for processor ',
253 & me,' set ',i2set(me)
257 call card_concat(weightcard)
258 call reada(weightcard,'WLONG',wlong,1.0D0)
259 call reada(weightcard,'WSC',wsc,wlong)
260 call reada(weightcard,'WSCP',wscp,wlong)
261 call reada(weightcard,'WELEC',welec,1.0D0)
262 call reada(weightcard,'WVDWPP',wvdwpp,welec)
263 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
264 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
265 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
266 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
267 call reada(weightcard,'WTURN3',wturn3,1.0D0)
268 call reada(weightcard,'WTURN4',wturn4,1.0D0)
269 call reada(weightcard,'WTURN6',wturn6,1.0D0)
270 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
271 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
272 call reada(weightcard,'WBOND',wbond,1.0D0)
273 call reada(weightcard,'WTOR',wtor,1.0D0)
274 call reada(weightcard,'WTORD',wtor_d,1.0D0)
275 call reada(weightcard,'WANG',wang,1.0D0)
276 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
277 call reada(weightcard,'SCAL14',scal14,0.4D0)
278 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
279 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
280 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
281 call reada(weightcard,'TEMP0',temp0,300.0d0)
282 if (index(weightcard,'SOFT').gt.0) ipot=6
283 C 12/1/95 Added weight for the multi-body term WCORR
284 call reada(weightcard,'WCORRH',wcorr,1.0D0)
285 if (wcorr4.gt.0.0d0) wcorr=wcorr4
293 hweights(i,7)=wel_loc
296 hweights(i,10)=wturn6
298 hweights(i,12)=wscloc
300 hweights(i,14)=wtor_d
301 hweights(i,15)=wstrain
302 hweights(i,16)=wvdwpp
304 hweights(i,18)=scal14
305 hweights(i,21)=wsccor
310 weights(i)=hweights(i2set(me),i)
334 call card_concat(weightcard)
335 call reada(weightcard,'WLONG',wlong,1.0D0)
336 call reada(weightcard,'WSC',wsc,wlong)
337 call reada(weightcard,'WSCP',wscp,wlong)
338 call reada(weightcard,'WELEC',welec,1.0D0)
339 call reada(weightcard,'WVDWPP',wvdwpp,welec)
340 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
341 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
342 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
343 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
344 call reada(weightcard,'WTURN3',wturn3,1.0D0)
345 call reada(weightcard,'WTURN4',wturn4,1.0D0)
346 call reada(weightcard,'WTURN6',wturn6,1.0D0)
347 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
348 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
349 call reada(weightcard,'WBOND',wbond,1.0D0)
350 call reada(weightcard,'WTOR',wtor,1.0D0)
351 call reada(weightcard,'WTORD',wtor_d,1.0D0)
352 call reada(weightcard,'WANG',wang,1.0D0)
353 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
354 call reada(weightcard,'SCAL14',scal14,0.4D0)
355 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
356 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
357 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
358 call reada(weightcard,'TEMP0',temp0,300.0d0)
359 if (index(weightcard,'SOFT').gt.0) ipot=6
360 C 12/1/95 Added weight for the multi-body term WCORR
361 call reada(weightcard,'WCORRH',wcorr,1.0D0)
362 if (wcorr4.gt.0.0d0) wcorr=wcorr4
384 if(me.eq.king.or..not.out1file)
385 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
386 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
388 10 format (/'Energy-term weights (unscaled):'//
389 & 'WSCC= ',f10.6,' (SC-SC)'/
390 & 'WSCP= ',f10.6,' (SC-p)'/
391 & 'WELEC= ',f10.6,' (p-p electr)'/
392 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
393 & 'WBOND= ',f10.6,' (stretching)'/
394 & 'WANG= ',f10.6,' (bending)'/
395 & 'WSCLOC= ',f10.6,' (SC local)'/
396 & 'WTOR= ',f10.6,' (torsional)'/
397 & 'WTORD= ',f10.6,' (double torsional)'/
398 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
399 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
400 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
401 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
402 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
403 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
404 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
405 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
406 & 'WTURN6= ',f10.6,' (turns, 6th order)')
407 if(me.eq.king.or..not.out1file)then
408 if (wcorr4.gt.0.0d0) then
409 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
410 & 'between contact pairs of peptide groups'
411 write (iout,'(2(a,f5.3/))')
412 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
413 & 'Range of quenching the correlation terms:',2*delt_corr
414 else if (wcorr.gt.0.0d0) then
415 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
416 & 'between contact pairs of peptide groups'
418 write (iout,'(a,f8.3)')
419 & 'Scaling factor of 1,4 SC-p interactions:',scal14
420 write (iout,'(a,f8.3)')
421 & 'General scaling factor of SC-p interactions:',scalscp
423 r0_corr=cutoff_corr-delt_corr
425 aad(i,1)=scalscp*aad(i,1)
426 aad(i,2)=scalscp*aad(i,2)
427 bad(i,1)=scalscp*bad(i,1)
428 bad(i,2)=scalscp*bad(i,2)
430 call rescale_weights(t_bath)
431 if(me.eq.king.or..not.out1file)
432 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
433 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
435 22 format (/'Energy-term weights (scaled):'//
436 & 'WSCC= ',f10.6,' (SC-SC)'/
437 & 'WSCP= ',f10.6,' (SC-p)'/
438 & 'WELEC= ',f10.6,' (p-p electr)'/
439 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
440 & 'WBOND= ',f10.6,' (stretching)'/
441 & 'WANG= ',f10.6,' (bending)'/
442 & 'WSCLOC= ',f10.6,' (SC local)'/
443 & 'WTOR= ',f10.6,' (torsional)'/
444 & 'WTORD= ',f10.6,' (double torsional)'/
445 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
446 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
447 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
448 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
449 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
450 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
451 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
452 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
453 & 'WTURN6= ',f10.6,' (turns, 6th order)')
454 if(me.eq.king.or..not.out1file)
455 & write (iout,*) "Reference temperature for weights calculation:",
457 call reada(weightcard,"D0CM",d0cm,3.78d0)
458 call reada(weightcard,"AKCM",akcm,15.1d0)
459 call reada(weightcard,"AKTH",akth,11.0d0)
460 call reada(weightcard,"AKCT",akct,12.0d0)
461 call reada(weightcard,"V1SS",v1ss,-1.08d0)
462 call reada(weightcard,"V2SS",v2ss,7.61d0)
463 call reada(weightcard,"V3SS",v3ss,13.7d0)
464 call reada(weightcard,"EBR",ebr,-5.50D0)
465 if(me.eq.king.or..not.out1file) then
466 write (iout,*) "Parameters of the SS-bond potential:"
467 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
469 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
470 write (iout,*) "EBR",ebr
471 print *,'indpdb=',indpdb,' pdbref=',pdbref
473 if (indpdb.gt.0 .or. pdbref) then
474 read(inp,'(a)') pdbname(1)
476 call contact(.false.,ncont_ref,icont_ref,co)
478 if (indpdb.eq.0) then
479 C Read sequence if not taken from the pdb file.
481 c print *,'nres=',nres
482 if (iscode.gt.0) then
483 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
485 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
487 C Convert sequence to numeric code
489 itype(i)=rescode(i,sequence(i),iscode)
491 C Assign initial virtual bond lengths
497 vbld(i+nres)=dsc(itype(i))
498 vbld_inv(i+nres)=dsc_inv(itype(i))
499 c write (iout,*) "i",i," itype",itype(i),
500 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
504 c print '(20i4)',(itype(i),i=1,nres)
507 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
509 if (itype(i).eq.21) then
513 else if (itype(i+1).ne.20) then
515 else if (itype(i).ne.20) then
522 if(me.eq.king.or..not.out1file)then
523 write (iout,*) "ITEL"
525 write (iout,*) i,itype(i),itel(i)
527 print *,'Call Read_Bridge.'
530 C 8/13/98 Set limits to generating the dihedral angles
535 read (inp,*) ndih_constr
536 if (ndih_constr.gt.0) then
538 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
539 if(me.eq.king.or..not.out1file)then
541 & 'There are',ndih_constr,' constraints on phi angles.'
543 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
547 phi0(i)=deg2rad*phi0(i)
548 drange(i)=deg2rad*drange(i)
550 if(me.eq.king.or..not.out1file)
551 & write (iout,*) 'FTORS',ftors
554 phibound(1,ii) = phi0(i)-drange(i)
555 phibound(2,ii) = phi0(i)+drange(i)
562 write (iout,'(a)') 'Boundaries in phi angle sampling:'
564 write (iout,'(a3,i5,2f10.1)')
565 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
571 cd print *,'NNT=',NNT,' NCT=',NCT
572 if (itype(1).eq.21) nnt=2
573 if (itype(nres).eq.21) nct=nct-1
575 if(me.eq.king.or..not.out1file)
576 & write (iout,'(a,i3)') 'nsup=',nsup
578 if (nsup.le.(nct-nnt+1)) then
579 do i=0,nct-nnt+1-nsup
580 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
586 & 'Error - sequences to be superposed do not match.'
589 do i=0,nsup-(nct-nnt+1)
590 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
592 nstart_sup=nstart_sup+i
598 & 'Error - sequences to be superposed do not match.'
601 if (nsup.eq.0) nsup=nct-nnt
602 if (nstart_sup.eq.0) nstart_sup=nnt
603 if (nstart_seq.eq.0) nstart_seq=nnt
604 if(me.eq.king.or..not.out1file)
605 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
606 & ' nstart_seq=',nstart_seq
608 c--- Zscore rms -------
609 if (nz_start.eq.0) nz_start=nnt
610 if (nz_end.eq.0 .and. nsup.gt.0) then
612 else if (nz_end.eq.0) then
615 if(me.eq.king.or..not.out1file)then
616 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
617 write (iout,*) 'IZ_SC=',iz_sc
619 c----------------------
622 if (.not.pdbref) then
623 call read_angles(inp,*38)
625 38 write (iout,'(a)') 'Error reading reference structure.'
627 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
628 stop 'Error reading reference structure'
632 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
641 call contact(.true.,ncont_ref,icont_ref,co)
643 if(me.eq.king.or..not.out1file)
644 & write (iout,*) 'Contact order:',co
646 if(me.eq.king.or..not.out1file)
647 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
650 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
652 if(me.eq.king.or..not.out1file)
653 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
654 & icont_ref(1,i),' ',
655 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
659 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
660 if (constr_dist.gt.0) then
661 call read_dist_constr
663 if (nhpb.gt.0) call hpb_partition
664 c write (iout,*) "After read_dist_constr nhpb",nhpb
666 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
667 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
668 & modecalc.ne.10) then
669 C If input structure hasn't been supplied from the PDB file read or generate
671 if (iranconf.eq.0 .and. .not. extconf) then
672 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
673 & write (iout,'(a)') 'Initial geometry will be read in.'
675 read(inp,'(8f10.5)',end=36,err=36)
676 & ((c(l,k),l=1,3),k=1,nres),
677 & ((c(l,k+nres),l=1,3),k=nnt,nct)
678 call int_from_cart1(.false.)
681 dc(j,i)=c(j,i+1)-c(j,i)
682 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
686 if (itype(i).ne.10) then
688 dc(j,i+nres)=c(j,i+nres)-c(j,i)
689 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
695 call read_angles(inp,*36)
698 36 write (iout,'(a)') 'Error reading angle file.'
700 call mpi_finalize( MPI_COMM_WORLD,IERR )
702 stop 'Error reading angle file.'
704 else if (extconf) then
705 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
706 & write (iout,'(a)') 'Extended chain initial geometry.'
708 theta(i)=90d0*deg2rad
714 alph(i)=110d0*deg2rad
717 omeg(i)=-120d0*deg2rad
720 if(me.eq.king.or..not.out1file)
721 & write (iout,'(a)') 'Random-generated initial geometry.'
725 if (me.eq.king .or. fg_rank.eq.0 .and. (
726 & modecalc.eq.12 .or. modecalc.eq.14) ) then
730 call gen_rand_conf(itmp,*30)
732 30 write (iout,*) 'Failed to generate random conformation',
734 write (*,*) 'Processor:',me,
735 & ' Failed to generate random conformation',
745 write (iout,'(a,i3,a)') 'Processor:',me,
746 & ' error in generating random conformation.'
747 write (*,'(a,i3,a)') 'Processor:',me,
748 & ' error in generating random conformation.'
751 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
758 elseif (modecalc.eq.4) then
759 read (inp,'(a)') intinname
760 open (intin,file=intinname,status='old',err=333)
761 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
762 & write (iout,'(a)') 'intinname',intinname
763 write (*,'(a)') 'Processor',myrank,' intinname',intinname
765 333 write (iout,'(2a)') 'Error opening angle file ',intinname
767 call MPI_Finalize(MPI_COMM_WORLD,IERR)
769 stop 'Error opening angle file.'
773 C Generate distance constraints, if the PDB structure is to be regularized.
774 if (nthread.gt.0) then
778 if (me.eq.king .or. .not. out1file)
780 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
781 write (iout,'(/a,i3,a)')
782 & 'The chain contains',ns,' disulfide-bridging cysteines.'
783 write (iout,'(20i4)') (iss(i),i=1,ns)
784 write (iout,'(/a/)') 'Pre-formed links are:'
790 if (me.eq.king.or..not.out1file)
791 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
792 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
797 if (i2ndstr.gt.0) call secstrp2dihc
798 c call geom_to_var(nvar,x)
799 c call etotal(energia(0))
800 c call enerprint(energia(0))
801 c call briefout(0,etot)
803 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
804 cd write (iout,'(a)') 'Variable list:'
805 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
806 if (indpdb.gt.0) then
808 read(1,'(a)',end=555,err=555) pdbname(i)
811 write (iout,*) "PDB files"
813 write (iout,'(i6,a)') i,pdbname(i)(:ilen(pdbname(i)))
817 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
818 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
819 & 'Processor',myrank,': end reading molecular data.'
823 c--------------------------------------------------------------------------
824 logical function seq_comp(itypea,itypeb,length)
826 integer length,itypea(length),itypeb(length)
829 if (itypea(i).ne.itypeb(i)) then
837 c-----------------------------------------------------------------------------
838 subroutine read_bridge
839 C Read information about disulfide bridges.
840 implicit real*8 (a-h,o-z)
845 include 'COMMON.IOUNITS'
848 include 'COMMON.INTERACT'
849 include 'COMMON.LOCAL'
850 include 'COMMON.NAMES'
851 include 'COMMON.CHAIN'
852 include 'COMMON.FFIELD'
853 include 'COMMON.SBRIDGE'
854 include 'COMMON.HEADER'
855 include 'COMMON.CONTROL'
856 include 'COMMON.DBASE'
857 include 'COMMON.THREAD'
858 include 'COMMON.TIME1'
859 include 'COMMON.SETUP'
860 C Read bridging residues.
861 read (inp,*) ns,(iss(i),i=1,ns)
863 if(me.eq.king.or..not.out1file)
864 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
865 C Check whether the specified bridging residues are cystines.
867 if (itype(iss(i)).ne.1) then
868 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
869 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
870 & ' can form a disulfide bridge?!!!'
871 write (*,'(2a,i3,a)')
872 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
873 & ' can form a disulfide bridge?!!!'
875 call MPI_Finalize(MPI_COMM_WORLD,ierror)
880 C Read preformed bridges.
882 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
883 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
886 C Check if the residues involved in bridges are in the specified list of
890 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
891 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
892 write (iout,'(a,i3,a)') 'Disulfide pair',i,
893 & ' contains residues present in other pairs.'
894 write (*,'(a,i3,a)') 'Disulfide pair',i,
895 & ' contains residues present in other pairs.'
897 call MPI_Finalize(MPI_COMM_WORLD,ierror)
903 if (ihpb(i).eq.iss(j)) goto 10
905 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
908 if (jhpb(i).eq.iss(j)) goto 20
910 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
923 c----------------------------------------------------------------------------
924 subroutine read_x(kanal,*)
925 implicit real*8 (a-h,o-z)
929 include 'COMMON.CHAIN'
930 include 'COMMON.IOUNITS'
931 include 'COMMON.CONTROL'
932 include 'COMMON.LOCAL'
933 include 'COMMON.INTERACT'
934 c Read coordinates from input
936 read(kanal,'(8f10.5)',end=10,err=10)
937 & ((c(l,k),l=1,3),k=1,nres),
938 & ((c(l,k+nres),l=1,3),k=nnt,nct)
941 c(j,2*nres)=c(j,nres)
943 call int_from_cart1(.false.)
946 dc(j,i)=c(j,i+1)-c(j,i)
947 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
951 if (itype(i).ne.10) then
953 dc(j,i+nres)=c(j,i+nres)-c(j,i)
954 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
962 c----------------------------------------------------------------------------
963 subroutine read_threadbase
964 implicit real*8 (a-h,o-z)
966 include 'COMMON.IOUNITS'
969 include 'COMMON.INTERACT'
970 include 'COMMON.LOCAL'
971 include 'COMMON.NAMES'
972 include 'COMMON.CHAIN'
973 include 'COMMON.FFIELD'
974 include 'COMMON.SBRIDGE'
975 include 'COMMON.HEADER'
976 include 'COMMON.CONTROL'
977 include 'COMMON.DBASE'
978 include 'COMMON.THREAD'
979 include 'COMMON.TIME1'
980 C Read pattern database for threading.
983 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
984 & nres_base(2,i),nres_base(3,i)
985 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
987 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
988 c & nres_base(2,i),nres_base(3,i)
989 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
993 if (weidis.eq.0.0D0) weidis=0.1D0
1002 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1003 write (iout,'(a,i5)') 'nexcl: ',nexcl
1004 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1007 c------------------------------------------------------------------------------
1008 subroutine setup_var
1009 implicit real*8 (a-h,o-z)
1010 include 'DIMENSIONS'
1011 include 'COMMON.IOUNITS'
1012 include 'COMMON.GEO'
1013 include 'COMMON.VAR'
1014 include 'COMMON.INTERACT'
1015 include 'COMMON.LOCAL'
1016 include 'COMMON.NAMES'
1017 include 'COMMON.CHAIN'
1018 include 'COMMON.FFIELD'
1019 include 'COMMON.SBRIDGE'
1020 include 'COMMON.HEADER'
1021 include 'COMMON.CONTROL'
1022 include 'COMMON.DBASE'
1023 include 'COMMON.THREAD'
1024 include 'COMMON.TIME1'
1025 C Set up variable list.
1031 if (itype(i).ne.10) then
1033 ialph(i,1)=nvar+nside
1037 if (indphi.gt.0) then
1039 else if (indback.gt.0) then
1044 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1047 c----------------------------------------------------------------------------
1048 subroutine gen_dist_constr
1049 C Generate CA distance constraints.
1050 implicit real*8 (a-h,o-z)
1051 include 'DIMENSIONS'
1052 include 'COMMON.IOUNITS'
1053 include 'COMMON.GEO'
1054 include 'COMMON.VAR'
1055 include 'COMMON.INTERACT'
1056 include 'COMMON.LOCAL'
1057 include 'COMMON.NAMES'
1058 include 'COMMON.CHAIN'
1059 include 'COMMON.FFIELD'
1060 include 'COMMON.SBRIDGE'
1061 include 'COMMON.HEADER'
1062 include 'COMMON.CONTROL'
1063 include 'COMMON.DBASE'
1064 include 'COMMON.THREAD'
1065 include 'COMMON.TIME1'
1066 dimension itype_pdb(maxres)
1067 common /pizda/ itype_pdb
1069 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1070 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1071 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1073 do i=nstart_sup,nstart_sup+nsup-1
1074 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1075 cd & ' seq_pdb', restyp(itype_pdb(i))
1076 do j=i+2,nstart_sup+nsup-1
1078 ihpb(nhpb)=i+nstart_seq-nstart_sup
1079 jhpb(nhpb)=j+nstart_seq-nstart_sup
1081 dhpb(nhpb)=dist(i,j)
1084 cd write (iout,'(a)') 'Distance constraints:'
1089 cd if (ii.gt.nres) then
1094 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1095 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1096 cd & dhpb(i),forcon(i)
1100 c----------------------------------------------------------------------------
1102 implicit real*8 (a-h,o-z)
1103 include 'DIMENSIONS'
1104 include 'COMMON.MAP'
1105 include 'COMMON.IOUNITS'
1106 character*3 angid(4) /'THE','PHI','ALP','OME'/
1107 character*80 mapcard,ucase
1109 read (inp,'(a)') mapcard
1110 mapcard=ucase(mapcard)
1111 if (index(mapcard,'PHI').gt.0) then
1113 else if (index(mapcard,'THE').gt.0) then
1115 else if (index(mapcard,'ALP').gt.0) then
1117 else if (index(mapcard,'OME').gt.0) then
1120 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1121 stop 'Error - illegal variable spec in MAP card.'
1123 call readi (mapcard,'RES1',res1(imap),0)
1124 call readi (mapcard,'RES2',res2(imap),0)
1125 if (res1(imap).eq.0) then
1126 res1(imap)=res2(imap)
1127 else if (res2(imap).eq.0) then
1128 res2(imap)=res1(imap)
1130 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1132 & 'Error - illegal definition of variable group in MAP.'
1133 stop 'Error - illegal definition of variable group in MAP.'
1135 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1136 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1137 call readi(mapcard,'NSTEP',nstep(imap),0)
1138 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1140 & 'Illegal boundary and/or step size specification in MAP.'
1141 stop 'Illegal boundary and/or step size specification in MAP.'
1146 c----------------------------------------------------------------------------
1147 subroutine read_minim
1148 implicit real*8 (a-h,o-z)
1149 include 'DIMENSIONS'
1150 include 'COMMON.MINIM'
1151 include 'COMMON.IOUNITS'
1153 character*320 minimcard
1154 call card_concat(minimcard)
1155 call readi(minimcard,'MAXMIN',maxmin,2000)
1156 call readi(minimcard,'MAXFUN',maxfun,5000)
1157 call readi(minimcard,'MINMIN',minmin,maxmin)
1158 call readi(minimcard,'MINFUN',minfun,maxmin)
1159 call reada(minimcard,'TOLF',tolf,1.0D-2)
1160 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1161 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1162 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1163 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1164 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1165 & 'Options in energy minimization:'
1166 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1167 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1168 & 'MinMin:',MinMin,' MinFun:',MinFun,
1169 & ' TolF:',TolF,' RTolF:',RTolF
1172 c----------------------------------------------------------------------------
1173 subroutine read_angles(kanal,*)
1174 implicit real*8 (a-h,o-z)
1175 include 'DIMENSIONS'
1176 include 'COMMON.GEO'
1177 include 'COMMON.VAR'
1178 include 'COMMON.CHAIN'
1179 include 'COMMON.IOUNITS'
1180 include 'COMMON.CONTROL'
1181 c Read angles from input
1183 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1184 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1185 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1186 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1189 c 9/7/01 avoid 180 deg valence angle
1190 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1192 theta(i)=deg2rad*theta(i)
1193 phi(i)=deg2rad*phi(i)
1194 alph(i)=deg2rad*alph(i)
1195 omeg(i)=deg2rad*omeg(i)
1200 c----------------------------------------------------------------------------
1201 subroutine reada(rekord,lancuch,wartosc,default)
1203 character*(*) rekord,lancuch
1204 double precision wartosc,default
1207 iread=index(rekord,lancuch)
1208 if (iread.eq.0) then
1212 iread=iread+ilen(lancuch)+1
1213 read (rekord(iread:),*,err=10,end=10) wartosc
1218 c----------------------------------------------------------------------------
1219 subroutine readi(rekord,lancuch,wartosc,default)
1221 character*(*) rekord,lancuch
1222 integer wartosc,default
1225 iread=index(rekord,lancuch)
1226 if (iread.eq.0) then
1230 iread=iread+ilen(lancuch)+1
1231 read (rekord(iread:),*,err=10,end=10) wartosc
1236 c----------------------------------------------------------------------------
1237 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1240 integer tablica(dim),default
1241 character*(*) rekord,lancuch
1248 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1249 if (iread.eq.0) return
1250 iread=iread+ilen(lancuch)+1
1251 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1254 c----------------------------------------------------------------------------
1255 subroutine multreada(rekord,lancuch,tablica,dim,default)
1258 double precision tablica(dim),default
1259 character*(*) rekord,lancuch
1266 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1267 if (iread.eq.0) return
1268 iread=iread+ilen(lancuch)+1
1269 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1272 c----------------------------------------------------------------------------
1273 subroutine card_concat(card)
1274 implicit real*8 (a-h,o-z)
1275 include 'DIMENSIONS'
1276 include 'COMMON.IOUNITS'
1278 character*80 karta,ucase
1280 read (inp,'(a)') karta
1283 do while (karta(80:80).eq.'&')
1284 card=card(:ilen(card)+1)//karta(:79)
1285 read (inp,'(a)') karta
1288 card=card(:ilen(card)+1)//karta
1291 c---------------------------------------------------------------------------------
1292 subroutine read_fragments
1293 implicit real*8 (a-h,o-z)
1294 include 'DIMENSIONS'
1298 include 'COMMON.SETUP'
1299 include 'COMMON.CHAIN'
1300 include 'COMMON.IOUNITS'
1302 include 'COMMON.CONTROL'
1303 read(inp,*) nset,nfrag,npair,nfrag_back
1304 if(me.eq.king.or..not.out1file)
1305 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
1306 & " nfrag_back",nfrag_back
1308 read(inp,*) mset(iset)
1310 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
1312 if(me.eq.king.or..not.out1file)
1313 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
1314 & ifrag(2,i,iset), qinfrag(i,iset)
1317 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
1319 if(me.eq.king.or..not.out1file)
1320 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
1321 & ipair(2,i,iset), qinpair(i,iset)
1324 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
1325 & wfrag_back(3,i,iset),
1326 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
1327 if(me.eq.king.or..not.out1file)
1328 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
1329 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
1334 c-------------------------------------------------------------------------------
1335 subroutine read_dist_constr
1336 implicit real*8 (a-h,o-z)
1337 include 'DIMENSIONS'
1341 include 'COMMON.SETUP'
1342 include 'COMMON.CONTROL'
1343 include 'COMMON.CHAIN'
1344 include 'COMMON.IOUNITS'
1345 include 'COMMON.SBRIDGE'
1346 integer ifrag_(2,100),ipair_(2,100)
1347 double precision wfrag_(100),wpair_(100)
1348 character*500 controlcard
1349 c write (iout,*) "Calling read_dist_constr"
1350 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1352 call card_concat(controlcard)
1353 call readi(controlcard,"NFRAG",nfrag_,0)
1354 call readi(controlcard,"NPAIR",npair_,0)
1355 call readi(controlcard,"NDIST",ndist_,0)
1356 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1357 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1358 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1359 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1360 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1361 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1362 c write (iout,*) "IFRAG"
1364 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1366 c write (iout,*) "IPAIR"
1368 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1370 if (.not.refstr .and. nfrag.gt.0) then
1372 & "ERROR: no reference structure to compute distance restraints"
1374 & "Restraints must be specified explicitly (NDIST=number)"
1377 if (nfrag.lt.2 .and. npair.gt.0) then
1378 write (iout,*) "ERROR: Less than 2 fragments specified",
1379 & " but distance restraints between pairs requested"
1384 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1385 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1386 & ifrag_(2,i)=nstart_sup+nsup-1
1387 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1389 if (wfrag_(i).gt.0.0d0) then
1390 do j=ifrag_(1,i),ifrag_(2,i)-1
1391 do k=j+1,ifrag_(2,i)
1392 write (iout,*) "j",j," k",k
1394 if (constr_dist.eq.1) then
1399 forcon(nhpb)=wfrag_(i)
1400 else if (constr_dist.eq.2) then
1401 if (ddjk.le.dist_cut) then
1406 forcon(nhpb)=wfrag_(i)
1413 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1416 if (.not.out1file .or. me.eq.king)
1417 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1418 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1420 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1421 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1428 if (wpair_(i).gt.0.0d0) then
1436 do j=ifrag_(1,ii),ifrag_(2,ii)
1437 do k=ifrag_(1,jj),ifrag_(2,jj)
1441 forcon(nhpb)=wpair_(i)
1442 dhpb(nhpb)=dist(j,k)
1444 if (.not.out1file .or. me.eq.king)
1445 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1446 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1448 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1449 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1456 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
1457 & ibecarb(i),forcon(nhpb+1)
1458 if (forcon(nhpb+1).gt.0.0d0) then
1460 if (ibecarb(i).gt.0) then
1461 ihpb(i)=ihpb(i)+nres
1462 jhpb(i)=jhpb(i)+nres
1464 if (dhpb(nhpb).eq.0.0d0)
1465 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1469 if (.not.out1file .or. me.eq.king) then
1472 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
1473 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
1481 c-------------------------------------------------------------------------------
1483 subroutine flush(iu)
1488 subroutine flush(iu)
1493 c------------------------------------------------------------------------------
1494 subroutine copy_to_tmp(source)
1495 include "DIMENSIONS"
1496 include "COMMON.IOUNITS"
1497 character*(*) source
1498 character* 256 tmpfile
1502 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1503 inquire(file=tmpfile,exist=ex)
1505 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1506 & " to temporary directory..."
1507 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1508 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1512 c------------------------------------------------------------------------------
1513 subroutine move_from_tmp(source)
1514 include "DIMENSIONS"
1515 include "COMMON.IOUNITS"
1516 character*(*) source
1519 write (*,*) "Moving ",source(:ilen(source)),
1520 & " from temporary directory to working directory"
1521 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1522 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1525 c------------------------------------------------------------------------------
1526 subroutine random_init(seed)
1528 C Initialize random number generator
1530 implicit real*8 (a-h,o-z)
1531 include 'DIMENSIONS'
1537 logical OKRandom, prng_restart
1539 integer iseed_array(4)
1541 include 'COMMON.IOUNITS'
1542 include 'COMMON.TIME1'
1543 include 'COMMON.THREAD'
1544 include 'COMMON.SBRIDGE'
1545 include 'COMMON.CONTROL'
1546 include 'COMMON.MCM'
1547 include 'COMMON.MAP'
1548 include 'COMMON.HEADER'
1549 csa include 'COMMON.CSA'
1550 include 'COMMON.CHAIN'
1551 include 'COMMON.MUCA'
1553 include 'COMMON.FFIELD'
1554 include 'COMMON.SETUP'
1555 iseed=-dint(dabs(seed))
1556 if (iseed.eq.0) then
1557 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1558 & 'Random seed undefined. The program will stop.'
1559 write (*,'(/80(1h*)/20x,a/80(1h*))')
1560 & 'Random seed undefined. The program will stop.'
1562 call mpi_finalize(mpi_comm_world,ierr)
1564 stop 'Bad random seed.'
1567 if (fg_rank.eq.0) then
1571 if(me.eq.king .or. .not. out1file)
1572 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1573 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1574 OKRandom = prng_restart(me,iseedi8)
1577 tmp=65536.0d0**(4-i)
1578 iseed_array(i) = dint(seed/tmp)
1579 seed=seed-iseed_array(i)*tmp
1581 if(me.eq.king .or. .not. out1file)
1582 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1583 & (iseed_array(i),i=1,4)
1584 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1585 & (iseed_array(i),i=1,4)
1586 OKRandom = prng_restart(me,iseed_array)
1589 r1=ran_number(0.0D0,1.0D0)
1590 if(me.eq.king .or. .not. out1file)
1591 & write (iout,*) 'ran_num',r1
1592 if (r1.lt.0.0d0) OKRandom=.false.
1594 if (.not.OKRandom) then
1595 write (iout,*) 'PRNG IS NOT WORKING!!!'
1596 print *,'PRNG IS NOT WORKING!!!'
1599 call mpi_abort(mpi_comm_world,error_msg,ierr)
1602 write (iout,*) 'too many processors for parallel prng'
1603 write (*,*) 'too many processors for parallel prng'
1611 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1615 c--------------------------------------------------------------------------
1616 subroutine pdbinput(iconf)
1617 implicit real*8 (a-h,o-z)
1618 include 'DIMENSIONS'
1623 include 'COMMON.IOUNITS'
1624 include 'COMMON.GEO'
1625 include 'COMMON.VAR'
1626 include 'COMMON.INTERACT'
1627 include 'COMMON.LOCAL'
1628 include 'COMMON.NAMES'
1629 include 'COMMON.CHAIN'
1630 include 'COMMON.FFIELD'
1631 include 'COMMON.SBRIDGE'
1632 include 'COMMON.HEADER'
1633 include 'COMMON.CONTROL'
1634 include 'COMMON.DBASE'
1635 include 'COMMON.THREAD'
1636 include 'COMMON.CONTACTS'
1637 include 'COMMON.TORCNSTR'
1638 include 'COMMON.TIME1'
1639 include 'COMMON.BOUNDS'
1641 include 'COMMON.REMD'
1642 include 'COMMON.SETUP'
1643 character*4 sequence(maxres)
1645 double precision x(maxvar)
1646 dimension itype_pdb(maxres)
1647 common /pizda/ itype_pdb
1648 logical seq_comp,fail
1649 double precision energia(0:n_ene)
1652 if(me.eq.king.or..not.out1file)
1653 & write (iout,'(2a)') 'PDB data will be read from file ',
1654 & pdbname(iconf)(:ilen(pdbname(iconf)))
1655 open(ipdbin,file=pdbname(iconf),status='old',err=33)
1657 33 write (iout,'(a)') 'Error opening PDB file.'
1660 c print *,'Begin reading pdb data'
1662 c print *,'Finished reading pdb data'
1663 if(me.eq.king.or..not.out1file)
1664 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
1665 & ' nstart_sup=',nstart_sup
1667 itype_pdb(i)=itype(i)
1671 nct=nstart_sup+nsup-1
1674 C Following 2 lines for diagnostics; comment out if not needed
1675 write (iout,*) "Before sideadd"
1677 if (me.eq.king.or..not.out1file)
1678 & write(iout,*)'Adding sidechains'
1685 do while (fail.and.nsi.le.maxsi)
1686 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
1689 if(fail) write(iout,*)'Adding sidechain failed for res ',
1690 & i,' after ',nsi,' trials'
1693 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
1695 C Following 2 lines for diagnostics; comment out if not needed
1696 write (iout,*) "After sideadd"