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 MCM control parameters if required
19 c if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 c if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 c if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 c if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 c if (modecalc.eq.8) then
32 c inquire (file="fort.40",exist=file_exist)
33 c if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 include 'COMMON.SBRIDGE'
70 include 'COMMON.CONTROL'
72 include 'COMMON.HEADER'
73 include 'COMMON.CHAIN'
74 include 'COMMON.FFIELD'
75 include 'COMMON.SETUP'
76 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
77 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
79 character*320 controlcard
84 read (INP,'(a)') titel
85 call card_concat(controlcard)
86 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
87 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
88 call reada(controlcard,'SEED',seed,0.0D0)
89 call random_init(seed)
90 C Set up the time limit (caution! The time must be input in minutes!)
91 read_cart=index(controlcard,'READ_CART').gt.0
92 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
93 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
94 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
95 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
96 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
97 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
98 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
99 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
100 call reada(controlcard,'DRMS',drms,0.1D0)
101 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
102 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
103 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
104 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
105 write (iout,'(a,f10.1)')'DRMS = ',drms
106 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
107 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
109 call readi(controlcard,'NZ_START',nz_start,0)
110 call readi(controlcard,'NZ_END',nz_end,0)
111 call readi(controlcard,'IZ_SC',iz_sc,0)
113 safety = 60.0d0*safety
116 call reada(controlcard,"T_BATH",t_bath,300.0d0)
117 minim=(index(controlcard,'MINIMIZE').gt.0)
118 dccart=(index(controlcard,'CART').gt.0)
119 overlapsc=(index(controlcard,'OVERLAP').gt.0)
120 overlapsc=.not.overlapsc
121 searchsc=(index(controlcard,'SEARCHSC').gt.0)
122 sideadd=(index(controlcard,'SIDEADD').gt.0)
123 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
124 outpdb=(index(controlcard,'PDBOUT').gt.0)
125 outmol2=(index(controlcard,'MOL2OUT').gt.0)
126 pdbref=(index(controlcard,'PDBREF').gt.0)
127 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
128 indpdb=index(controlcard,'PDBSTART')
129 extconf=(index(controlcard,'EXTCONF').gt.0)
130 call readi(controlcard,'IPRINT',iprint,0)
131 call readi(controlcard,'MAXGEN',maxgen,10000)
132 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
133 call readi(controlcard,"KDIAG",kdiag,0)
134 call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
135 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
136 & write (iout,*) "RESCALE_MODE",rescale_mode
137 split_ene=index(controlcard,'SPLIT_ENE').gt.0
138 if (index(controlcard,'REGULAR').gt.0.0D0) then
139 call reada(controlcard,'WEIDIS',weidis,0.1D0)
143 if (index(controlcard,'CHECKGRAD').gt.0) then
145 if (index(controlcard,'CART').gt.0) then
147 elseif (index(controlcard,'CARINT').gt.0) then
152 elseif (index(controlcard,'THREAD').gt.0) then
154 call readi(controlcard,'THREAD',nthread,0)
155 if (nthread.gt.0) then
156 call reada(controlcard,'WEIDIS',weidis,0.1D0)
159 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
160 stop 'Error termination in Read_Control.'
162 else if (index(controlcard,'MCMA').gt.0) then
164 else if (index(controlcard,'MCEE').gt.0) then
166 else if (index(controlcard,'MULTCONF').gt.0) then
168 else if (index(controlcard,'MAP').gt.0) then
170 call readi(controlcard,'MAP',nmap,0)
171 else if (index(controlcard,'CSA').gt.0) then
173 crc else if (index(controlcard,'ZSCORE').gt.0) then
175 crc ZSCORE is rm from UNRES, modecalc=9 is available
178 cfcm else if (index(controlcard,'MCMF').gt.0) then
180 else if (index(controlcard,'SOFTREG').gt.0) then
182 else if (index(controlcard,'CHECK_BOND').gt.0) then
184 else if (index(controlcard,'TEST').gt.0) then
186 else if (index(controlcard,'MD').gt.0) then
188 else if (index(controlcard,'RE ').gt.0) then
192 lmuca=index(controlcard,'MUCA').gt.0
193 call readi(controlcard,'MUCADYN',mucadyn,0)
194 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
195 if (lmuca .and. (me.eq.king .or. .not.out1file ))
197 write (iout,*) 'MUCADYN=',mucadyn
198 write (iout,*) 'MUCASMOOTH=',muca_smooth
201 iscode=index(controlcard,'ONE_LETTER')
202 indphi=index(controlcard,'PHI')
203 indback=index(controlcard,'BACK')
204 iranconf=index(controlcard,'RAND_CONF')
205 i2ndstr=index(controlcard,'USE_SEC_PRED')
206 gradout=index(controlcard,'GRADOUT').gt.0
207 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
209 if(me.eq.king.or..not.out1file)
210 & write (iout,'(2a)') diagmeth(kdiag),
211 & ' routine used to diagonalize matrices.'
214 c------------------------------------------------------------------------------
217 C Read molecular data.
219 implicit real*8 (a-h,o-z)
225 include 'COMMON.IOUNITS'
228 include 'COMMON.INTERACT'
229 include 'COMMON.LOCAL'
230 include 'COMMON.NAMES'
231 include 'COMMON.CHAIN'
232 include 'COMMON.FFIELD'
233 include 'COMMON.SBRIDGE'
234 include 'COMMON.HEADER'
235 include 'COMMON.CONTROL'
236 c include 'COMMON.DBASE'
237 c include 'COMMON.THREAD'
238 include 'COMMON.CONTACTS'
239 include 'COMMON.TORCNSTR'
240 include 'COMMON.TIME1'
241 include 'COMMON.BOUNDS'
242 c include 'COMMON.MD'
243 c include 'COMMON.REMD'
244 include 'COMMON.SETUP'
245 character*4 sequence(maxres)
247 double precision x(maxvar)
248 character*256 pdbfile
249 character*320 weightcard
250 character*80 weightcard_t,ucase
251 dimension itype_pdb(maxres)
252 common /pizda/ itype_pdb
253 logical seq_comp,fail
254 double precision energia(0:n_ene)
260 C Read weights of the subsequent energy terms.
262 call card_concat(weightcard)
263 call reada(weightcard,'WLONG',wlong,1.0D0)
264 call reada(weightcard,'WSC',wsc,wlong)
265 call reada(weightcard,'WSCP',wscp,wlong)
266 call reada(weightcard,'WELEC',welec,1.0D0)
267 call reada(weightcard,'WVDWPP',wvdwpp,welec)
268 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
269 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
270 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
271 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
272 call reada(weightcard,'WTURN3',wturn3,1.0D0)
273 call reada(weightcard,'WTURN4',wturn4,1.0D0)
274 call reada(weightcard,'WTURN6',wturn6,1.0D0)
275 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
276 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
277 call reada(weightcard,'WBOND',wbond,1.0D0)
278 call reada(weightcard,'WTOR',wtor,1.0D0)
279 call reada(weightcard,'WTORD',wtor_d,1.0D0)
280 call reada(weightcard,'WANG',wang,1.0D0)
281 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
282 call reada(weightcard,'SCAL14',scal14,0.4D0)
283 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
284 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
285 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
286 call reada(weightcard,'TEMP0',temp0,300.0d0)
287 if (index(weightcard,'SOFT').gt.0) ipot=6
288 C 12/1/95 Added weight for the multi-body term WCORR
289 call reada(weightcard,'WCORRH',wcorr,1.0D0)
290 if (wcorr4.gt.0.0d0) wcorr=wcorr4
311 if(me.eq.king.or..not.out1file)
312 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
313 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
315 10 format (/'Energy-term weights (unscaled):'//
316 & 'WSCC= ',f10.6,' (SC-SC)'/
317 & 'WSCP= ',f10.6,' (SC-p)'/
318 & 'WELEC= ',f10.6,' (p-p electr)'/
319 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
320 & 'WBOND= ',f10.6,' (stretching)'/
321 & 'WANG= ',f10.6,' (bending)'/
322 & 'WSCLOC= ',f10.6,' (SC local)'/
323 & 'WTOR= ',f10.6,' (torsional)'/
324 & 'WTORD= ',f10.6,' (double torsional)'/
325 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
326 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
327 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
328 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
329 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
330 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
331 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
332 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
333 & 'WTURN6= ',f10.6,' (turns, 6th order)')
334 if(me.eq.king.or..not.out1file)then
335 if (wcorr4.gt.0.0d0) then
336 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
337 & 'between contact pairs of peptide groups'
338 write (iout,'(2(a,f5.3/))')
339 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
340 & 'Range of quenching the correlation terms:',2*delt_corr
341 else if (wcorr.gt.0.0d0) then
342 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
343 & 'between contact pairs of peptide groups'
345 write (iout,'(a,f8.3)')
346 & 'Scaling factor of 1,4 SC-p interactions:',scal14
347 write (iout,'(a,f8.3)')
348 & 'General scaling factor of SC-p interactions:',scalscp
350 r0_corr=cutoff_corr-delt_corr
352 aad(i,1)=scalscp*aad(i,1)
353 aad(i,2)=scalscp*aad(i,2)
354 bad(i,1)=scalscp*bad(i,1)
355 bad(i,2)=scalscp*bad(i,2)
357 c call rescale_weights(t_bath)
358 if(me.eq.king.or..not.out1file)
359 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
360 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
362 22 format (/'Energy-term weights (scaled):'//
363 & 'WSCC= ',f10.6,' (SC-SC)'/
364 & 'WSCP= ',f10.6,' (SC-p)'/
365 & 'WELEC= ',f10.6,' (p-p electr)'/
366 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
367 & 'WBOND= ',f10.6,' (stretching)'/
368 & 'WANG= ',f10.6,' (bending)'/
369 & 'WSCLOC= ',f10.6,' (SC local)'/
370 & 'WTOR= ',f10.6,' (torsional)'/
371 & 'WTORD= ',f10.6,' (double torsional)'/
372 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
373 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
374 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
375 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
376 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
377 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
378 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
379 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
380 & 'WTURN6= ',f10.6,' (turns, 6th order)')
381 if(me.eq.king.or..not.out1file)
382 & write (iout,*) "Reference temperature for weights calculation:",
384 call reada(weightcard,"D0CM",d0cm,3.78d0)
385 call reada(weightcard,"AKCM",akcm,15.1d0)
386 call reada(weightcard,"AKTH",akth,11.0d0)
387 call reada(weightcard,"AKCT",akct,12.0d0)
388 call reada(weightcard,"V1SS",v1ss,-1.08d0)
389 call reada(weightcard,"V2SS",v2ss,7.61d0)
390 call reada(weightcard,"V3SS",v3ss,13.7d0)
391 call reada(weightcard,"EBR",ebr,-5.50D0)
392 if(me.eq.king.or..not.out1file) then
393 write (iout,*) "Parameters of the SS-bond potential:"
394 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
396 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
397 write (iout,*) "EBR",ebr
398 print *,'indpdb=',indpdb,' pdbref=',pdbref
400 if (indpdb.gt.0 .or. pdbref) then
401 read(inp,'(a)') pdbfile
402 if(me.eq.king.or..not.out1file)
403 & write (iout,'(2a)') 'PDB data will be read from file ',
404 & pdbfile(:ilen(pdbfile))
405 open(ipdbin,file=pdbfile,status='old',err=33)
407 33 write (iout,'(a)') 'Error opening PDB file.'
410 c print *,'Begin reading pdb data'
412 c print *,'Finished reading pdb data'
413 if(me.eq.king.or..not.out1file)
414 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
415 & ' nstart_sup=',nstart_sup
417 itype_pdb(i)=itype(i)
421 nct=nstart_sup+nsup-1
422 c call contact(.false.,ncont_ref,icont_ref,co)
425 if(me.eq.king.or..not.out1file)
426 & write(iout,*)'Adding sidechains'
433 do while (fail.and.nsi.le.maxsi)
434 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
437 if(fail) write(iout,*)'Adding sidechain failed for res ',
438 & i,' after ',nsi,' trials'
443 if (indpdb.eq.0) then
444 C Read sequence if not taken from the pdb file.
446 c print *,'nres=',nres
447 if (iscode.gt.0) then
448 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
450 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
452 C Convert sequence to numeric code
454 itype(i)=rescode(i,sequence(i),iscode)
456 C Assign initial virtual bond lengths
462 vbld(i+nres)=dsc(itype(i))
463 vbld_inv(i+nres)=dsc_inv(itype(i))
464 c write (iout,*) "i",i," itype",itype(i),
465 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
469 c print '(20i4)',(itype(i),i=1,nres)
472 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
474 if (itype(i).eq.21) then
478 else if (itype(i+1).ne.20) then
480 else if (itype(i).ne.20) then
487 if(me.eq.king.or..not.out1file)then
488 write (iout,*) "ITEL"
490 write (iout,*) i,itype(i),itel(i)
492 print *,'Call Read_Bridge.'
495 C 8/13/98 Set limits to generating the dihedral angles
500 read (inp,*) ndih_constr
501 if (ndih_constr.gt.0) then
503 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
504 if(me.eq.king.or..not.out1file)then
506 & 'There are',ndih_constr,' constraints on phi angles.'
508 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
512 phi0(i)=deg2rad*phi0(i)
513 drange(i)=deg2rad*drange(i)
515 if(me.eq.king.or..not.out1file)
516 & write (iout,*) 'FTORS',ftors
519 phibound(1,ii) = phi0(i)-drange(i)
520 phibound(2,ii) = phi0(i)+drange(i)
527 write (iout,'(a)') 'Boundaries in phi angle sampling:'
529 write (iout,'(a3,i5,2f10.1)')
530 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
536 cd print *,'NNT=',NNT,' NCT=',NCT
537 if (itype(1).eq.21) nnt=2
538 if (itype(nres).eq.21) nct=nct-1
540 if(me.eq.king.or..not.out1file)
541 & write (iout,'(a,i3)') 'nsup=',nsup
543 if (nsup.le.(nct-nnt+1)) then
544 do i=0,nct-nnt+1-nsup
545 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
551 & 'Error - sequences to be superposed do not match.'
554 do i=0,nsup-(nct-nnt+1)
555 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
557 nstart_sup=nstart_sup+i
563 & 'Error - sequences to be superposed do not match.'
566 if (nsup.eq.0) nsup=nct-nnt
567 if (nstart_sup.eq.0) nstart_sup=nnt
568 if (nstart_seq.eq.0) nstart_seq=nnt
569 if(me.eq.king.or..not.out1file)
570 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
571 & ' nstart_seq=',nstart_seq
573 c--- Zscore rms -------
574 if (nz_start.eq.0) nz_start=nnt
575 if (nz_end.eq.0 .and. nsup.gt.0) then
577 else if (nz_end.eq.0) then
580 if(me.eq.king.or..not.out1file)then
581 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
582 write (iout,*) 'IZ_SC=',iz_sc
584 c----------------------
587 if (.not.pdbref) then
588 call read_angles(inp,*38)
590 38 write (iout,'(a)') 'Error reading reference structure.'
592 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
593 stop 'Error reading reference structure'
597 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
606 c call contact(.true.,ncont_ref,icont_ref,co)
608 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
610 if (constr_dist.gt.0) call read_dist_constr
611 c write (iout,*) "After read_dist_constr nhpb",nhpb
613 if(me.eq.king.or..not.out1file)
614 & write (iout,*) 'Contact order:',co
616 if(me.eq.king.or..not.out1file)
617 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
620 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
622 if(me.eq.king.or..not.out1file)
623 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
624 & icont_ref(1,i),' ',
625 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
629 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
630 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
631 & modecalc.ne.10) then
632 C If input structure hasn't been supplied from the PDB file read or generate
634 if (iranconf.eq.0 .and. .not. extconf) then
635 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
636 & write (iout,'(a)') 'Initial geometry will be read in.'
638 read(inp,'(8f10.5)',end=36,err=36)
639 & ((c(l,k),l=1,3),k=1,nres),
640 & ((c(l,k+nres),l=1,3),k=nnt,nct)
641 call int_from_cart1(.false.)
644 dc(j,i)=c(j,i+1)-c(j,i)
645 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
649 if (itype(i).ne.10) then
651 dc(j,i+nres)=c(j,i+nres)-c(j,i)
652 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
658 call read_angles(inp,*36)
661 36 write (iout,'(a)') 'Error reading angle file.'
663 call mpi_finalize( MPI_COMM_WORLD,IERR )
665 stop 'Error reading angle file.'
667 else if (extconf) then
668 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
669 & write (iout,'(a)') 'Extended chain initial geometry.'
671 theta(i)=90d0*deg2rad
677 alph(i)=110d0*deg2rad
680 omeg(i)=-120d0*deg2rad
683 if(me.eq.king.or..not.out1file)
684 & write (iout,'(a)') 'Random-generated initial geometry.'
688 if (me.eq.king .or. fg_rank.eq.0 .and. (
689 & modecalc.eq.12 .or. modecalc.eq.14) ) then
693 call gen_rand_conf(itmp,*30)
695 30 write (iout,*) 'Failed to generate random conformation',
697 write (*,*) 'Processor:',me,
698 & ' Failed to generate random conformation',
708 write (iout,'(a,i3,a)') 'Processor:',me,
709 & ' error in generating random conformation.'
710 write (*,'(a,i3,a)') 'Processor:',me,
711 & ' error in generating random conformation.'
714 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
720 call gen_rand_conf(itmp,*31)
722 31 write (iout,*) 'Failed to generate random conformation',
724 write (*,*) 'Failed to generate random conformation',
727 write (iout,'(a,i3,a)') 'Processor:',me,
728 & ' error in generating random conformation.'
729 write (*,'(a,i3,a)') 'Processor:',me,
730 & ' error in generating random conformation.'
735 elseif (modecalc.eq.4) then
736 read (inp,'(a)') intinname
737 open (intin,file=intinname,status='old',err=333)
738 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
739 & write (iout,'(a)') 'intinname',intinname
740 write (*,'(a)') 'Processor',myrank,' intinname',intinname
742 333 write (iout,'(2a)') 'Error opening angle file ',intinname
744 call MPI_Finalize(MPI_COMM_WORLD,IERR)
746 stop 'Error opening angle file.'
750 C Generate distance constraints, if the PDB structure is to be regularized.
751 c if (nthread.gt.0) then
752 c call read_threadbase
755 if (me.eq.king .or. .not. out1file)
757 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
758 write (iout,'(/a,i3,a)')
759 & 'The chain contains',ns,' disulfide-bridging cysteines.'
760 write (iout,'(20i4)') (iss(i),i=1,ns)
761 write (iout,'(/a/)') 'Pre-formed links are:'
767 if (me.eq.king.or..not.out1file)
768 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
769 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
774 c if (i2ndstr.gt.0) call secstrp2dihc
775 c call geom_to_var(nvar,x)
776 c call etotal(energia(0))
777 c call enerprint(energia(0))
778 c call briefout(0,etot)
780 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
781 cd write (iout,'(a)') 'Variable list:'
782 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
784 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
785 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
786 & 'Processor',myrank,': end reading molecular data.'
790 c--------------------------------------------------------------------------
791 logical function seq_comp(itypea,itypeb,length)
793 integer length,itypea(length),itypeb(length)
796 if (itypea(i).ne.itypeb(i)) then
804 c-----------------------------------------------------------------------------
805 subroutine read_bridge
806 C Read information about disulfide bridges.
807 implicit real*8 (a-h,o-z)
812 include 'COMMON.IOUNITS'
815 include 'COMMON.INTERACT'
816 include 'COMMON.LOCAL'
817 include 'COMMON.NAMES'
818 include 'COMMON.CHAIN'
819 include 'COMMON.FFIELD'
820 include 'COMMON.SBRIDGE'
821 include 'COMMON.HEADER'
822 include 'COMMON.CONTROL'
823 c include 'COMMON.DBASE'
824 c include 'COMMON.THREAD'
825 include 'COMMON.TIME1'
826 include 'COMMON.SETUP'
827 C Read bridging residues.
828 read (inp,*) ns,(iss(i),i=1,ns)
830 if(me.eq.king.or..not.out1file)
831 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
832 C Check whether the specified bridging residues are cystines.
834 if (itype(iss(i)).ne.1) then
835 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
836 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
837 & ' can form a disulfide bridge?!!!'
838 write (*,'(2a,i3,a)')
839 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
840 & ' can form a disulfide bridge?!!!'
842 call MPI_Finalize(MPI_COMM_WORLD,ierror)
847 C Read preformed bridges.
849 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
850 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
853 C Check if the residues involved in bridges are in the specified list of
857 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
858 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
859 write (iout,'(a,i3,a)') 'Disulfide pair',i,
860 & ' contains residues present in other pairs.'
861 write (*,'(a,i3,a)') 'Disulfide pair',i,
862 & ' contains residues present in other pairs.'
864 call MPI_Finalize(MPI_COMM_WORLD,ierror)
870 if (ihpb(i).eq.iss(j)) goto 10
872 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
875 if (jhpb(i).eq.iss(j)) goto 20
877 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
890 c----------------------------------------------------------------------------
891 subroutine read_x(kanal,*)
892 implicit real*8 (a-h,o-z)
896 include 'COMMON.CHAIN'
897 include 'COMMON.IOUNITS'
898 include 'COMMON.CONTROL'
899 include 'COMMON.LOCAL'
900 include 'COMMON.INTERACT'
901 c Read coordinates from input
903 read(kanal,'(8f10.5)',end=10,err=10)
904 & ((c(l,k),l=1,3),k=1,nres),
905 & ((c(l,k+nres),l=1,3),k=nnt,nct)
908 c(j,2*nres)=c(j,nres)
910 call int_from_cart1(.false.)
913 dc(j,i)=c(j,i+1)-c(j,i)
914 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
918 if (itype(i).ne.10) then
920 dc(j,i+nres)=c(j,i+nres)-c(j,i)
921 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
929 c------------------------------------------------------------------------------
931 implicit real*8 (a-h,o-z)
933 include 'COMMON.IOUNITS'
936 include 'COMMON.INTERACT'
937 include 'COMMON.LOCAL'
938 include 'COMMON.NAMES'
939 include 'COMMON.CHAIN'
940 include 'COMMON.FFIELD'
941 include 'COMMON.SBRIDGE'
942 include 'COMMON.HEADER'
943 include 'COMMON.CONTROL'
944 c include 'COMMON.DBASE'
945 c include 'COMMON.THREAD'
946 include 'COMMON.TIME1'
947 C Set up variable list.
953 if (itype(i).ne.10) then
955 ialph(i,1)=nvar+nside
959 if (indphi.gt.0) then
961 else if (indback.gt.0) then
966 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
969 c----------------------------------------------------------------------------
970 subroutine gen_dist_constr
971 C Generate CA distance constraints.
972 implicit real*8 (a-h,o-z)
974 include 'COMMON.IOUNITS'
977 include 'COMMON.INTERACT'
978 include 'COMMON.LOCAL'
979 include 'COMMON.NAMES'
980 include 'COMMON.CHAIN'
981 include 'COMMON.FFIELD'
982 include 'COMMON.SBRIDGE'
983 include 'COMMON.HEADER'
984 include 'COMMON.CONTROL'
985 c include 'COMMON.DBASE'
986 c include 'COMMON.THREAD'
987 include 'COMMON.TIME1'
988 dimension itype_pdb(maxres)
989 common /pizda/ itype_pdb
991 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
992 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
993 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
995 do i=nstart_sup,nstart_sup+nsup-1
996 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
997 cd & ' seq_pdb', restyp(itype_pdb(i))
998 do j=i+2,nstart_sup+nsup-1
1000 ihpb(nhpb)=i+nstart_seq-nstart_sup
1001 jhpb(nhpb)=j+nstart_seq-nstart_sup
1003 dhpb(nhpb)=dist(i,j)
1006 cd write (iout,'(a)') 'Distance constraints:'
1011 cd if (ii.gt.nres) then
1016 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1017 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1018 cd & dhpb(i),forcon(i)
1023 subroutine read_minim
1024 implicit real*8 (a-h,o-z)
1025 include 'DIMENSIONS'
1026 include 'COMMON.MINIM'
1027 include 'COMMON.IOUNITS'
1029 character*320 minimcard
1030 call card_concat(minimcard)
1031 call readi(minimcard,'MAXMIN',maxmin,2000)
1032 call readi(minimcard,'MAXFUN',maxfun,5000)
1033 call readi(minimcard,'MINMIN',minmin,maxmin)
1034 call readi(minimcard,'MINFUN',minfun,maxmin)
1035 call reada(minimcard,'TOLF',tolf,1.0D-2)
1036 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1037 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1038 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1039 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1040 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1041 & 'Options in energy minimization:'
1042 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1043 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1044 & 'MinMin:',MinMin,' MinFun:',MinFun,
1045 & ' TolF:',TolF,' RTolF:',RTolF
1048 c----------------------------------------------------------------------------
1049 subroutine read_angles(kanal,*)
1050 implicit real*8 (a-h,o-z)
1051 include 'DIMENSIONS'
1052 include 'COMMON.GEO'
1053 include 'COMMON.VAR'
1054 include 'COMMON.CHAIN'
1055 include 'COMMON.IOUNITS'
1056 include 'COMMON.CONTROL'
1057 c Read angles from input
1059 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1060 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1061 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1062 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1065 c 9/7/01 avoid 180 deg valence angle
1066 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1068 theta(i)=deg2rad*theta(i)
1069 phi(i)=deg2rad*phi(i)
1070 alph(i)=deg2rad*alph(i)
1071 omeg(i)=deg2rad*omeg(i)
1076 c----------------------------------------------------------------------------
1077 subroutine reada(rekord,lancuch,wartosc,default)
1079 character*(*) rekord,lancuch
1080 double precision wartosc,default
1083 iread=index(rekord,lancuch)
1084 if (iread.eq.0) then
1088 iread=iread+ilen(lancuch)+1
1089 read (rekord(iread:),*,err=10,end=10) wartosc
1094 c----------------------------------------------------------------------------
1095 subroutine readi(rekord,lancuch,wartosc,default)
1097 character*(*) rekord,lancuch
1098 integer wartosc,default
1101 iread=index(rekord,lancuch)
1102 if (iread.eq.0) then
1106 iread=iread+ilen(lancuch)+1
1107 read (rekord(iread:),*,err=10,end=10) wartosc
1112 c----------------------------------------------------------------------------
1113 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1116 integer tablica(dim),default
1117 character*(*) rekord,lancuch
1124 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1125 if (iread.eq.0) return
1126 iread=iread+ilen(lancuch)+1
1127 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1130 c----------------------------------------------------------------------------
1131 subroutine multreada(rekord,lancuch,tablica,dim,default)
1134 double precision tablica(dim),default
1135 character*(*) rekord,lancuch
1142 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1143 if (iread.eq.0) return
1144 iread=iread+ilen(lancuch)+1
1145 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1148 c----------------------------------------------------------------------------
1149 subroutine openunits
1150 implicit real*8 (a-h,o-z)
1151 include 'DIMENSIONS'
1154 character*16 form,nodename
1157 include 'COMMON.SETUP'
1158 include 'COMMON.IOUNITS'
1159 c include 'COMMON.MD'
1160 include 'COMMON.CONTROL'
1161 integer lenpre,lenpot,ilen,lentmp
1163 character*3 out1file_text,ucase
1166 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1167 call getenv_loc("PREFIX",prefix)
1169 call getenv_loc("POT",pot)
1170 call getenv_loc("DIRTMP",tmpdir)
1171 call getenv_loc("CURDIR",curdir)
1172 call getenv_loc("OUT1FILE",out1file_text)
1173 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1174 out1file_text=ucase(out1file_text)
1175 if (out1file_text(1:1).eq."Y") then
1178 out1file=fg_rank.gt.0
1183 if (lentmp.gt.0) then
1184 write (*,'(80(1h!))')
1185 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1186 write (*,'(80(1h!))')
1187 write (*,*)"All output files will be on node /tmp directory."
1189 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1190 if (me.eq.king) then
1191 write (*,*) "The master node is ",nodename
1192 else if (fg_rank.eq.0) then
1193 write (*,*) "I am the CG slave node ",nodename
1195 write (*,*) "I am the FG slave node ",nodename
1198 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1199 lenpre = lentmp+lenpre+1
1201 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1202 C Get the names and open the input files
1203 #if defined(WINIFL) || defined(WINPGI)
1204 open(1,file=pref_orig(:ilen(pref_orig))//
1205 & '.inp',status='old',readonly,shared)
1206 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1207 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1208 C Get parameter filenames and open the parameter files.
1209 call getenv_loc('BONDPAR',bondname)
1210 open (ibond,file=bondname,status='old',readonly,shared)
1211 call getenv_loc('THETPAR',thetname)
1212 open (ithep,file=thetname,status='old',readonly,shared)
1214 call getenv_loc('THETPARPDB',thetname_pdb)
1215 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1217 call getenv_loc('ROTPAR',rotname)
1218 open (irotam,file=rotname,status='old',readonly,shared)
1220 call getenv_loc('ROTPARPDB',rotname_pdb)
1221 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1223 call getenv_loc('TORPAR',torname)
1224 open (itorp,file=torname,status='old',readonly,shared)
1225 call getenv_loc('TORDPAR',tordname)
1226 open (itordp,file=tordname,status='old',readonly,shared)
1227 call getenv_loc('FOURIER',fouriername)
1228 open (ifourier,file=fouriername,status='old',readonly,shared)
1229 call getenv_loc('ELEPAR',elename)
1230 open (ielep,file=elename,status='old',readonly,shared)
1231 call getenv_loc('SIDEPAR',sidename)
1232 open (isidep,file=sidename,status='old',readonly,shared)
1233 #elif (defined CRAY) || (defined AIX)
1234 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1236 c print *,"Processor",myrank," opened file 1"
1237 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1238 c print *,"Processor",myrank," opened file 9"
1239 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1240 C Get parameter filenames and open the parameter files.
1241 call getenv_loc('BONDPAR',bondname)
1242 open (ibond,file=bondname,status='old',action='read')
1243 c print *,"Processor",myrank," opened file IBOND"
1244 call getenv_loc('THETPAR',thetname)
1245 open (ithep,file=thetname,status='old',action='read')
1246 c print *,"Processor",myrank," opened file ITHEP"
1248 call getenv_loc('THETPARPDB',thetname_pdb)
1249 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1251 call getenv_loc('ROTPAR',rotname)
1252 open (irotam,file=rotname,status='old',action='read')
1253 c print *,"Processor",myrank," opened file IROTAM"
1255 call getenv_loc('ROTPARPDB',rotname_pdb)
1256 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1258 call getenv_loc('TORPAR',torname)
1259 open (itorp,file=torname,status='old',action='read')
1260 c print *,"Processor",myrank," opened file ITORP"
1261 call getenv_loc('TORDPAR',tordname)
1262 open (itordp,file=tordname,status='old',action='read')
1263 c print *,"Processor",myrank," opened file ITORDP"
1264 call getenv_loc('SCCORPAR',sccorname)
1265 open (isccor,file=sccorname,status='old',action='read')
1266 c print *,"Processor",myrank," opened file ISCCOR"
1267 call getenv_loc('FOURIER',fouriername)
1268 open (ifourier,file=fouriername,status='old',action='read')
1269 c print *,"Processor",myrank," opened file IFOURIER"
1270 call getenv_loc('ELEPAR',elename)
1271 open (ielep,file=elename,status='old',action='read')
1272 c print *,"Processor",myrank," opened file IELEP"
1273 call getenv_loc('SIDEPAR',sidename)
1274 open (isidep,file=sidename,status='old',action='read')
1275 c print *,"Processor",myrank," opened file ISIDEP"
1276 c print *,"Processor",myrank," opened parameter files"
1278 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1279 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1280 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1281 C Get parameter filenames and open the parameter files.
1282 call getenv_loc('BONDPAR',bondname)
1283 open (ibond,file=bondname,status='old')
1284 call getenv_loc('THETPAR',thetname)
1285 open (ithep,file=thetname,status='old')
1287 call getenv_loc('THETPARPDB',thetname_pdb)
1288 open (ithep_pdb,file=thetname_pdb,status='old')
1290 call getenv_loc('ROTPAR',rotname)
1291 open (irotam,file=rotname,status='old')
1293 call getenv_loc('ROTPARPDB',rotname_pdb)
1294 open (irotam_pdb,file=rotname_pdb,status='old')
1296 call getenv_loc('TORPAR',torname)
1297 open (itorp,file=torname,status='old')
1298 call getenv_loc('TORDPAR',tordname)
1299 open (itordp,file=tordname,status='old')
1300 call getenv_loc('SCCORPAR',sccorname)
1301 open (isccor,file=sccorname,status='old')
1302 call getenv_loc('FOURIER',fouriername)
1303 open (ifourier,file=fouriername,status='old')
1304 call getenv_loc('ELEPAR',elename)
1305 open (ielep,file=elename,status='old')
1306 call getenv_loc('SIDEPAR',sidename)
1307 open (isidep,file=sidename,status='old')
1309 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1311 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1312 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1313 C Get parameter filenames and open the parameter files.
1314 call getenv_loc('BONDPAR',bondname)
1315 open (ibond,file=bondname,status='old',readonly)
1316 call getenv_loc('THETPAR',thetname)
1317 open (ithep,file=thetname,status='old',readonly)
1319 call getenv_loc('THETPARPDB',thetname_pdb)
1320 print *,"thetname_pdb ",thetname_pdb
1321 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1322 print *,ithep_pdb," opened"
1324 call getenv_loc('ROTPAR',rotname)
1325 open (irotam,file=rotname,status='old',readonly)
1327 call getenv_loc('ROTPARPDB',rotname_pdb)
1328 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1330 call getenv_loc('TORPAR',torname)
1331 open (itorp,file=torname,status='old',readonly)
1332 call getenv_loc('TORDPAR',tordname)
1333 open (itordp,file=tordname,status='old',readonly)
1334 call getenv_loc('SCCORPAR',sccorname)
1335 open (isccor,file=sccorname,status='old',readonly)
1336 call getenv_loc('FOURIER',fouriername)
1337 open (ifourier,file=fouriername,status='old',readonly)
1338 call getenv_loc('ELEPAR',elename)
1339 open (ielep,file=elename,status='old',readonly)
1340 call getenv_loc('SIDEPAR',sidename)
1341 open (isidep,file=sidename,status='old',readonly)
1345 C 8/9/01 In the newest version SCp interaction constants are read from a file
1346 C Use -DOLDSCP to use hard-coded constants instead.
1348 call getenv_loc('SCPPAR',scpname)
1349 #if defined(WINIFL) || defined(WINPGI)
1350 open (iscpp,file=scpname,status='old',readonly,shared)
1351 #elif (defined CRAY) || (defined AIX)
1352 open (iscpp,file=scpname,status='old',action='read')
1354 open (iscpp,file=scpname,status='old')
1356 open (iscpp,file=scpname,status='old',readonly)
1359 call getenv_loc('PATTERN',patname)
1360 #if defined(WINIFL) || defined(WINPGI)
1361 open (icbase,file=patname,status='old',readonly,shared)
1362 #elif (defined CRAY) || (defined AIX)
1363 open (icbase,file=patname,status='old',action='read')
1365 open (icbase,file=patname,status='old')
1367 open (icbase,file=patname,status='old',readonly)
1370 C Open output file only for CG processes
1371 c print *,"Processor",myrank," fg_rank",fg_rank
1372 if (fg_rank.eq.0) then
1374 if (nodes.eq.1) then
1377 npos = dlog10(dfloat(nodes-1))+1
1379 if (npos.lt.3) npos=3
1380 write (liczba,'(i1)') npos
1381 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1383 write (liczba,form) me
1384 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1385 & liczba(:ilen(liczba))
1386 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1388 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1390 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1391 & liczba(:ilen(liczba))//'.mol2'
1392 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1393 & liczba(:ilen(liczba))//'.stat'
1395 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1396 & //liczba(:ilen(liczba))//'.stat')
1397 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1400 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1401 c & liczba(:ilen(liczba))//'.const'
1406 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1407 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1408 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1409 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1410 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1412 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1414 rest2name=prefix(:ilen(prefix))//'.rst'
1416 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1419 #if defined(AIX) || defined(PGI)
1420 if (me.eq.king .or. .not. out1file)
1421 & open(iout,file=outname,status='unknown')
1423 if (fg_rank.gt.0) then
1424 write (liczba,'(i3.3)') myrank/nfgtasks
1425 write (ll,'(bz,i3.3)') fg_rank
1426 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1431 open(igeom,file=intname,status='unknown',position='append')
1432 open(ipdb,file=pdbname,status='unknown')
1433 open(imol2,file=mol2name,status='unknown')
1434 open(istat,file=statname,status='unknown',position='append')
1436 c1out open(iout,file=outname,status='unknown')
1439 if (me.eq.king .or. .not.out1file)
1440 & open(iout,file=outname,status='unknown')
1442 if (fg_rank.gt.0) then
1443 write (liczba,'(i3.3)') myrank/nfgtasks
1444 write (ll,'(bz,i3.3)') fg_rank
1445 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1450 open(igeom,file=intname,status='unknown',access='append')
1451 open(ipdb,file=pdbname,status='unknown')
1452 open(imol2,file=mol2name,status='unknown')
1453 open(istat,file=statname,status='unknown',access='append')
1455 c1out open(iout,file=outname,status='unknown')
1458 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1459 csa_seed=prefix(:lenpre)//'.CSA.seed'
1460 csa_history=prefix(:lenpre)//'.CSA.history'
1461 csa_bank=prefix(:lenpre)//'.CSA.bank'
1462 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1463 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1464 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1465 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1466 csa_int=prefix(:lenpre)//'.int'
1467 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1468 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1469 csa_in=prefix(:lenpre)//'.CSA.in'
1470 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1473 write (iout,'(80(1h-))')
1474 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1475 write (iout,'(80(1h-))')
1476 write (iout,*) "Input file : ",
1477 & pref_orig(:ilen(pref_orig))//'.inp'
1478 write (iout,*) "Output file : ",
1479 & outname(:ilen(outname))
1481 write (iout,*) "Sidechain potential file : ",
1482 & sidename(:ilen(sidename))
1484 write (iout,*) "SCp potential file : ",
1485 & scpname(:ilen(scpname))
1487 write (iout,*) "Electrostatic potential file : ",
1488 & elename(:ilen(elename))
1489 write (iout,*) "Cumulant coefficient file : ",
1490 & fouriername(:ilen(fouriername))
1491 write (iout,*) "Torsional parameter file : ",
1492 & torname(:ilen(torname))
1493 write (iout,*) "Double torsional parameter file : ",
1494 & tordname(:ilen(tordname))
1495 write (iout,*) "SCCOR parameter file : ",
1496 & sccorname(:ilen(sccorname))
1497 write (iout,*) "Bond & inertia constant file : ",
1498 & bondname(:ilen(bondname))
1499 write (iout,*) "Bending parameter file : ",
1500 & thetname(:ilen(thetname))
1501 write (iout,*) "Rotamer parameter file : ",
1502 & rotname(:ilen(rotname))
1503 write (iout,*) "Threading database : ",
1504 & patname(:ilen(patname))
1506 &write (iout,*)" DIRTMP : ",
1508 write (iout,'(80(1h-))')
1512 c----------------------------------------------------------------------------
1513 subroutine card_concat(card)
1514 implicit real*8 (a-h,o-z)
1515 include 'DIMENSIONS'
1516 include 'COMMON.IOUNITS'
1518 character*80 karta,ucase
1520 read (inp,'(a)') karta
1523 do while (karta(80:80).eq.'&')
1524 card=card(:ilen(card)+1)//karta(:79)
1525 read (inp,'(a)') karta
1528 card=card(:ilen(card)+1)//karta
1532 subroutine read_dist_constr
1533 implicit real*8 (a-h,o-z)
1534 include 'DIMENSIONS'
1538 include 'COMMON.SETUP'
1539 include 'COMMON.CONTROL'
1540 include 'COMMON.CHAIN'
1541 include 'COMMON.IOUNITS'
1542 include 'COMMON.SBRIDGE'
1543 integer ifrag_(2,100),ipair_(2,100)
1544 double precision wfrag_(100),wpair_(100)
1545 character*500 controlcard
1546 c write (iout,*) "Calling read_dist_constr"
1547 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1549 call card_concat(controlcard)
1550 call readi(controlcard,"NFRAG",nfrag_,0)
1551 call readi(controlcard,"NPAIR",npair_,0)
1552 call readi(controlcard,"NDIST",ndist_,0)
1553 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1554 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1555 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1556 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1557 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1558 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1559 c write (iout,*) "IFRAG"
1561 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1563 c write (iout,*) "IPAIR"
1565 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1569 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1570 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1571 & ifrag_(2,i)=nstart_sup+nsup-1
1572 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1574 if (wfrag_(i).gt.0.0d0) then
1575 do j=ifrag_(1,i),ifrag_(2,i)-1
1576 do k=j+1,ifrag_(2,i)
1577 write (iout,*) "j",j," k",k
1579 if (constr_dist.eq.1) then
1584 forcon(nhpb)=wfrag_(i)
1585 else if (constr_dist.eq.2) then
1586 if (ddjk.le.dist_cut) then
1591 forcon(nhpb)=wfrag_(i)
1598 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1601 if (.not.out1file .or. me.eq.king)
1602 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1603 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1605 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1606 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1613 if (wpair_(i).gt.0.0d0) then
1621 do j=ifrag_(1,ii),ifrag_(2,ii)
1622 do k=ifrag_(1,jj),ifrag_(2,jj)
1626 forcon(nhpb)=wpair_(i)
1627 dhpb(nhpb)=dist(j,k)
1629 if (.not.out1file .or. me.eq.king)
1630 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1631 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1633 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1634 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1641 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1642 if (forcon(nhpb+1).gt.0.0d0) then
1644 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1646 if (.not.out1file .or. me.eq.king)
1647 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1648 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1650 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1651 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1658 c-------------------------------------------------------------------------------
1660 subroutine flush(iu)
1665 subroutine flush(iu)
1670 c------------------------------------------------------------------------------
1671 subroutine copy_to_tmp(source)
1672 include "DIMENSIONS"
1673 include "COMMON.IOUNITS"
1674 character*(*) source
1675 character* 256 tmpfile
1679 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1680 inquire(file=tmpfile,exist=ex)
1682 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1683 & " to temporary directory..."
1684 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1685 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1689 c------------------------------------------------------------------------------
1690 subroutine move_from_tmp(source)
1691 include "DIMENSIONS"
1692 include "COMMON.IOUNITS"
1693 character*(*) source
1696 write (*,*) "Moving ",source(:ilen(source)),
1697 & " from temporary directory to working directory"
1698 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1699 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1702 c------------------------------------------------------------------------------
1703 subroutine random_init(seed)
1705 C Initialize random number generator
1707 implicit real*8 (a-h,o-z)
1708 include 'DIMENSIONS'
1714 logical OKRandom, prng_restart
1716 integer iseed_array(4)
1718 include 'COMMON.IOUNITS'
1719 include 'COMMON.TIME1'
1720 c include 'COMMON.THREAD'
1721 include 'COMMON.SBRIDGE'
1722 include 'COMMON.CONTROL'
1723 include 'COMMON.MCM'
1724 c include 'COMMON.MAP'
1725 include 'COMMON.HEADER'
1726 c include 'COMMON.CSA'
1727 include 'COMMON.CHAIN'
1728 c include 'COMMON.MUCA'
1729 c include 'COMMON.MD'
1730 include 'COMMON.FFIELD'
1731 include 'COMMON.SETUP'
1732 iseed=-dint(dabs(seed))
1733 if (iseed.eq.0) then
1734 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1735 & 'Random seed undefined. The program will stop.'
1736 write (*,'(/80(1h*)/20x,a/80(1h*))')
1737 & 'Random seed undefined. The program will stop.'
1739 call mpi_finalize(mpi_comm_world,ierr)
1741 stop 'Bad random seed.'
1744 if (fg_rank.eq.0) then
1748 if(me.eq.king .or. .not. out1file)
1749 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1750 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1751 OKRandom = prng_restart(me,iseedi8)
1754 tmp=65536.0d0**(4-i)
1755 iseed_array(i) = dint(seed/tmp)
1756 seed=seed-iseed_array(i)*tmp
1758 if(me.eq.king .or. .not. out1file)
1759 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1760 & (iseed_array(i),i=1,4)
1761 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1762 & (iseed_array(i),i=1,4)
1763 OKRandom = prng_restart(me,iseed_array)
1766 c r1=ran_number(0.0D0,1.0D0)
1767 if(me.eq.king .or. .not. out1file)
1768 & write (iout,*) 'ran_num',r1
1769 if (r1.lt.0.0d0) OKRandom=.false.
1771 if (.not.OKRandom) then
1772 write (iout,*) 'PRNG IS NOT WORKING!!!'
1773 print *,'PRNG IS NOT WORKING!!!'
1776 call mpi_abort(mpi_comm_world,error_msg,ierr)
1779 write (iout,*) 'too many processors for parallel prng'
1780 write (*,*) 'too many processors for parallel prng'
1788 c write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)