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,'NOSEARCHSC').gt.0)
122 searchsc=.not.searchsc
123 sideadd=(index(controlcard,'SIDEADD').gt.0)
124 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
125 outpdb=(index(controlcard,'PDBOUT').gt.0)
126 outmol2=(index(controlcard,'MOL2OUT').gt.0)
127 pdbref=(index(controlcard,'PDBREF').gt.0)
128 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
129 indpdb=index(controlcard,'PDBSTART')
130 extconf=(index(controlcard,'EXTCONF').gt.0)
131 call readi(controlcard,'IPRINT',iprint,0)
132 call readi(controlcard,'MAXGEN',maxgen,10000)
133 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
134 call readi(controlcard,"KDIAG",kdiag,0)
135 call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
136 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
137 & write (iout,*) "RESCALE_MODE",rescale_mode
138 split_ene=index(controlcard,'SPLIT_ENE').gt.0
139 if (index(controlcard,'REGULAR').gt.0.0D0) then
140 call reada(controlcard,'WEIDIS',weidis,0.1D0)
144 if (index(controlcard,'CHECKGRAD').gt.0) then
146 if (index(controlcard,'CART').gt.0) then
148 elseif (index(controlcard,'CARINT').gt.0) then
153 elseif (index(controlcard,'THREAD').gt.0) then
155 call readi(controlcard,'THREAD',nthread,0)
156 if (nthread.gt.0) then
157 call reada(controlcard,'WEIDIS',weidis,0.1D0)
160 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
161 stop 'Error termination in Read_Control.'
163 else if (index(controlcard,'MCMA').gt.0) then
165 else if (index(controlcard,'MCEE').gt.0) then
167 else if (index(controlcard,'MULTCONF').gt.0) then
169 else if (index(controlcard,'MAP').gt.0) then
171 call readi(controlcard,'MAP',nmap,0)
172 else if (index(controlcard,'CSA').gt.0) then
174 crc else if (index(controlcard,'ZSCORE').gt.0) then
176 crc ZSCORE is rm from UNRES, modecalc=9 is available
179 cfcm else if (index(controlcard,'MCMF').gt.0) then
181 else if (index(controlcard,'SOFTREG').gt.0) then
183 else if (index(controlcard,'CHECK_BOND').gt.0) then
185 else if (index(controlcard,'TEST').gt.0) then
187 else if (index(controlcard,'MD').gt.0) then
189 else if (index(controlcard,'RE ').gt.0) then
193 lmuca=index(controlcard,'MUCA').gt.0
194 call readi(controlcard,'MUCADYN',mucadyn,0)
195 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
196 if (lmuca .and. (me.eq.king .or. .not.out1file ))
198 write (iout,*) 'MUCADYN=',mucadyn
199 write (iout,*) 'MUCASMOOTH=',muca_smooth
202 iscode=index(controlcard,'ONE_LETTER')
203 indphi=index(controlcard,'PHI')
204 indback=index(controlcard,'BACK')
205 iranconf=index(controlcard,'RAND_CONF')
206 i2ndstr=index(controlcard,'USE_SEC_PRED')
207 gradout=index(controlcard,'GRADOUT').gt.0
208 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
210 if(me.eq.king.or..not.out1file)
211 & write (iout,'(2a)') diagmeth(kdiag),
212 & ' routine used to diagonalize matrices.'
215 c------------------------------------------------------------------------------
218 C Read molecular data.
220 implicit real*8 (a-h,o-z)
226 include 'COMMON.IOUNITS'
229 include 'COMMON.INTERACT'
230 include 'COMMON.LOCAL'
231 include 'COMMON.NAMES'
232 include 'COMMON.CHAIN'
233 include 'COMMON.FFIELD'
234 include 'COMMON.SBRIDGE'
235 include 'COMMON.HEADER'
236 include 'COMMON.CONTROL'
237 c include 'COMMON.DBASE'
238 c include 'COMMON.THREAD'
239 include 'COMMON.CONTACTS'
240 include 'COMMON.TORCNSTR'
241 include 'COMMON.TIME1'
242 include 'COMMON.BOUNDS'
243 c include 'COMMON.MD'
244 c include 'COMMON.REMD'
245 include 'COMMON.SETUP'
246 character*4 sequence(maxres)
248 double precision x(maxvar)
249 character*256 pdbfile
250 character*320 weightcard
251 character*80 weightcard_t,ucase
252 dimension itype_pdb(maxres)
253 common /pizda/ itype_pdb
254 logical seq_comp,fail
255 double precision energia(0:n_ene)
261 C Read weights of the subsequent energy terms.
263 call card_concat(weightcard)
264 call reada(weightcard,'WLONG',wlong,1.0D0)
265 call reada(weightcard,'WSC',wsc,wlong)
266 call reada(weightcard,'WSCP',wscp,wlong)
267 call reada(weightcard,'WELEC',welec,1.0D0)
268 call reada(weightcard,'WVDWPP',wvdwpp,welec)
269 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
270 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
271 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
272 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
273 call reada(weightcard,'WTURN3',wturn3,1.0D0)
274 call reada(weightcard,'WTURN4',wturn4,1.0D0)
275 call reada(weightcard,'WTURN6',wturn6,1.0D0)
276 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
277 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
278 call reada(weightcard,'WBOND',wbond,1.0D0)
279 call reada(weightcard,'WTOR',wtor,1.0D0)
280 call reada(weightcard,'WTORD',wtor_d,1.0D0)
281 call reada(weightcard,'WANG',wang,1.0D0)
282 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
283 call reada(weightcard,'SCAL14',scal14,0.4D0)
284 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
285 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
286 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
287 call reada(weightcard,'TEMP0',temp0,300.0d0)
288 if (index(weightcard,'SOFT').gt.0) ipot=6
289 C 12/1/95 Added weight for the multi-body term WCORR
290 call reada(weightcard,'WCORRH',wcorr,1.0D0)
291 if (wcorr4.gt.0.0d0) wcorr=wcorr4
312 if(me.eq.king.or..not.out1file)
313 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
314 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
316 10 format (/'Energy-term weights (unscaled):'//
317 & 'WSCC= ',f10.6,' (SC-SC)'/
318 & 'WSCP= ',f10.6,' (SC-p)'/
319 & 'WELEC= ',f10.6,' (p-p electr)'/
320 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
321 & 'WBOND= ',f10.6,' (stretching)'/
322 & 'WANG= ',f10.6,' (bending)'/
323 & 'WSCLOC= ',f10.6,' (SC local)'/
324 & 'WTOR= ',f10.6,' (torsional)'/
325 & 'WTORD= ',f10.6,' (double torsional)'/
326 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
327 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
328 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
329 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
330 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
331 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
332 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
333 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
334 & 'WTURN6= ',f10.6,' (turns, 6th order)')
335 if(me.eq.king.or..not.out1file)then
336 if (wcorr4.gt.0.0d0) then
337 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
338 & 'between contact pairs of peptide groups'
339 write (iout,'(2(a,f5.3/))')
340 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
341 & 'Range of quenching the correlation terms:',2*delt_corr
342 else if (wcorr.gt.0.0d0) then
343 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
344 & 'between contact pairs of peptide groups'
346 write (iout,'(a,f8.3)')
347 & 'Scaling factor of 1,4 SC-p interactions:',scal14
348 write (iout,'(a,f8.3)')
349 & 'General scaling factor of SC-p interactions:',scalscp
351 r0_corr=cutoff_corr-delt_corr
353 aad(i,1)=scalscp*aad(i,1)
354 aad(i,2)=scalscp*aad(i,2)
355 bad(i,1)=scalscp*bad(i,1)
356 bad(i,2)=scalscp*bad(i,2)
358 c call rescale_weights(t_bath)
359 if(me.eq.king.or..not.out1file)
360 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
361 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
363 22 format (/'Energy-term weights (scaled):'//
364 & 'WSCC= ',f10.6,' (SC-SC)'/
365 & 'WSCP= ',f10.6,' (SC-p)'/
366 & 'WELEC= ',f10.6,' (p-p electr)'/
367 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
368 & 'WBOND= ',f10.6,' (stretching)'/
369 & 'WANG= ',f10.6,' (bending)'/
370 & 'WSCLOC= ',f10.6,' (SC local)'/
371 & 'WTOR= ',f10.6,' (torsional)'/
372 & 'WTORD= ',f10.6,' (double torsional)'/
373 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
374 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
375 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
376 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
377 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
378 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
379 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
380 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
381 & 'WTURN6= ',f10.6,' (turns, 6th order)')
382 if(me.eq.king.or..not.out1file)
383 & write (iout,*) "Reference temperature for weights calculation:",
385 call reada(weightcard,"D0CM",d0cm,3.78d0)
386 call reada(weightcard,"AKCM",akcm,15.1d0)
387 call reada(weightcard,"AKTH",akth,11.0d0)
388 call reada(weightcard,"AKCT",akct,12.0d0)
389 call reada(weightcard,"V1SS",v1ss,-1.08d0)
390 call reada(weightcard,"V2SS",v2ss,7.61d0)
391 call reada(weightcard,"V3SS",v3ss,13.7d0)
392 call reada(weightcard,"EBR",ebr,-5.50D0)
393 if(me.eq.king.or..not.out1file) then
394 write (iout,*) "Parameters of the SS-bond potential:"
395 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
397 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
398 write (iout,*) "EBR",ebr
399 print *,'indpdb=',indpdb,' pdbref=',pdbref
401 if (indpdb.gt.0 .or. pdbref) then
402 read(inp,'(a)') pdbfile
403 if(me.eq.king.or..not.out1file)
404 & write (iout,'(2a)') 'PDB data will be read from file ',
405 & pdbfile(:ilen(pdbfile))
406 open(ipdbin,file=pdbfile,status='old',err=33)
408 33 write (iout,'(a)') 'Error opening PDB file.'
411 c print *,'Begin reading pdb data'
413 c print *,'Finished reading pdb data'
414 if(me.eq.king.or..not.out1file)
415 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
416 & ' nstart_sup=',nstart_sup
418 itype_pdb(i)=itype(i)
422 nct=nstart_sup+nsup-1
423 c call contact(.false.,ncont_ref,icont_ref,co)
426 if(me.eq.king.or..not.out1file)
427 & write(iout,*)'Adding sidechains'
434 do while (fail.and.nsi.le.maxsi)
435 c call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
438 if(fail) write(iout,*)'Adding sidechain failed for res ',
439 & i,' after ',nsi,' trials'
444 if (indpdb.eq.0) then
445 C Read sequence if not taken from the pdb file.
447 c print *,'nres=',nres
448 if (iscode.gt.0) then
449 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
451 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
453 C Convert sequence to numeric code
455 itype(i)=rescode(i,sequence(i),iscode)
457 C Assign initial virtual bond lengths
463 vbld(i+nres)=dsc(itype(i))
464 vbld_inv(i+nres)=dsc_inv(itype(i))
465 c write (iout,*) "i",i," itype",itype(i),
466 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
470 c print '(20i4)',(itype(i),i=1,nres)
473 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
475 if (itype(i).eq.21) then
479 else if (itype(i+1).ne.20) then
481 else if (itype(i).ne.20) then
488 if(me.eq.king.or..not.out1file)then
489 write (iout,*) "ITEL"
491 write (iout,*) i,itype(i),itel(i)
493 print *,'Call Read_Bridge.'
496 C 8/13/98 Set limits to generating the dihedral angles
501 read (inp,*) ndih_constr
502 if (ndih_constr.gt.0) then
504 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
505 if(me.eq.king.or..not.out1file)then
507 & 'There are',ndih_constr,' constraints on phi angles.'
509 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
513 phi0(i)=deg2rad*phi0(i)
514 drange(i)=deg2rad*drange(i)
516 if(me.eq.king.or..not.out1file)
517 & write (iout,*) 'FTORS',ftors
520 phibound(1,ii) = phi0(i)-drange(i)
521 phibound(2,ii) = phi0(i)+drange(i)
528 write (iout,'(a)') 'Boundaries in phi angle sampling:'
530 write (iout,'(a3,i5,2f10.1)')
531 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
537 cd print *,'NNT=',NNT,' NCT=',NCT
538 if (itype(1).eq.21) nnt=2
539 if (itype(nres).eq.21) nct=nct-1
541 if(me.eq.king.or..not.out1file)
542 & write (iout,'(a,i3)') 'nsup=',nsup
544 if (nsup.le.(nct-nnt+1)) then
545 do i=0,nct-nnt+1-nsup
546 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
552 & 'Error - sequences to be superposed do not match.'
555 do i=0,nsup-(nct-nnt+1)
556 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
558 nstart_sup=nstart_sup+i
564 & 'Error - sequences to be superposed do not match.'
567 if (nsup.eq.0) nsup=nct-nnt
568 if (nstart_sup.eq.0) nstart_sup=nnt
569 if (nstart_seq.eq.0) nstart_seq=nnt
570 if(me.eq.king.or..not.out1file)
571 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
572 & ' nstart_seq=',nstart_seq
574 c--- Zscore rms -------
575 if (nz_start.eq.0) nz_start=nnt
576 if (nz_end.eq.0 .and. nsup.gt.0) then
578 else if (nz_end.eq.0) then
581 if(me.eq.king.or..not.out1file)then
582 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
583 write (iout,*) 'IZ_SC=',iz_sc
585 c----------------------
588 if (.not.pdbref) then
589 call read_angles(inp,*38)
591 38 write (iout,'(a)') 'Error reading reference structure.'
593 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
594 stop 'Error reading reference structure'
598 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
607 c call contact(.true.,ncont_ref,icont_ref,co)
609 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
611 if (constr_dist.gt.0) call read_dist_constr
612 c write (iout,*) "After read_dist_constr nhpb",nhpb
614 if(me.eq.king.or..not.out1file)
615 & write (iout,*) 'Contact order:',co
617 if(me.eq.king.or..not.out1file)
618 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
621 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
623 if(me.eq.king.or..not.out1file)
624 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
625 & icont_ref(1,i),' ',
626 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
630 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
631 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
632 & modecalc.ne.10) then
633 C If input structure hasn't been supplied from the PDB file read or generate
635 if (iranconf.eq.0 .and. .not. extconf) then
636 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
637 & write (iout,'(a)') 'Initial geometry will be read in.'
639 read(inp,'(8f10.5)',end=36,err=36)
640 & ((c(l,k),l=1,3),k=1,nres),
641 & ((c(l,k+nres),l=1,3),k=nnt,nct)
642 call int_from_cart1(.false.)
645 dc(j,i)=c(j,i+1)-c(j,i)
646 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
650 if (itype(i).ne.10) then
652 dc(j,i+nres)=c(j,i+nres)-c(j,i)
653 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
659 call read_angles(inp,*36)
662 36 write (iout,'(a)') 'Error reading angle file.'
664 call mpi_finalize( MPI_COMM_WORLD,IERR )
666 stop 'Error reading angle file.'
668 else if (extconf) then
669 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
670 & write (iout,'(a)') 'Extended chain initial geometry.'
672 theta(i)=90d0*deg2rad
678 alph(i)=110d0*deg2rad
681 omeg(i)=-120d0*deg2rad
684 if(me.eq.king.or..not.out1file)
685 & write (iout,'(a)') 'Random-generated initial geometry.'
689 if (me.eq.king .or. fg_rank.eq.0 .and. (
690 & modecalc.eq.12 .or. modecalc.eq.14) ) then
694 c call gen_rand_conf(itmp,*30)
696 30 write (iout,*) 'Failed to generate random conformation',
698 write (*,*) 'Processor:',me,
699 & ' Failed to generate random conformation',
709 write (iout,'(a,i3,a)') 'Processor:',me,
710 & ' error in generating random conformation.'
711 write (*,'(a,i3,a)') 'Processor:',me,
712 & ' error in generating random conformation.'
715 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
721 c call gen_rand_conf(itmp,*31)
723 31 write (iout,*) 'Failed to generate random conformation',
725 write (*,*) 'Failed to generate random conformation',
728 write (iout,'(a,i3,a)') 'Processor:',me,
729 & ' error in generating random conformation.'
730 write (*,'(a,i3,a)') 'Processor:',me,
731 & ' error in generating random conformation.'
736 elseif (modecalc.eq.4) then
737 read (inp,'(a)') intinname
738 open (intin,file=intinname,status='old',err=333)
739 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
740 & write (iout,'(a)') 'intinname',intinname
741 write (*,'(a)') 'Processor',myrank,' intinname',intinname
743 333 write (iout,'(2a)') 'Error opening angle file ',intinname
745 call MPI_Finalize(MPI_COMM_WORLD,IERR)
747 stop 'Error opening angle file.'
751 C Generate distance constraints, if the PDB structure is to be regularized.
752 c if (nthread.gt.0) then
753 c call read_threadbase
756 if (me.eq.king .or. .not. out1file)
758 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
759 write (iout,'(/a,i3,a)')
760 & 'The chain contains',ns,' disulfide-bridging cysteines.'
761 write (iout,'(20i4)') (iss(i),i=1,ns)
762 write (iout,'(/a/)') 'Pre-formed links are:'
768 if (me.eq.king.or..not.out1file)
769 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
770 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
775 c if (i2ndstr.gt.0) call secstrp2dihc
776 c call geom_to_var(nvar,x)
777 c call etotal(energia(0))
778 c call enerprint(energia(0))
779 c call briefout(0,etot)
781 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
782 cd write (iout,'(a)') 'Variable list:'
783 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
785 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
786 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
787 & 'Processor',myrank,': end reading molecular data.'
791 c--------------------------------------------------------------------------
792 logical function seq_comp(itypea,itypeb,length)
794 integer length,itypea(length),itypeb(length)
797 if (itypea(i).ne.itypeb(i)) then
805 c-----------------------------------------------------------------------------
806 subroutine read_bridge
807 C Read information about disulfide bridges.
808 implicit real*8 (a-h,o-z)
813 include 'COMMON.IOUNITS'
816 include 'COMMON.INTERACT'
817 include 'COMMON.LOCAL'
818 include 'COMMON.NAMES'
819 include 'COMMON.CHAIN'
820 include 'COMMON.FFIELD'
821 include 'COMMON.SBRIDGE'
822 include 'COMMON.HEADER'
823 include 'COMMON.CONTROL'
824 c include 'COMMON.DBASE'
825 c include 'COMMON.THREAD'
826 include 'COMMON.TIME1'
827 include 'COMMON.SETUP'
828 C Read bridging residues.
829 read (inp,*) ns,(iss(i),i=1,ns)
831 if(me.eq.king.or..not.out1file)
832 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
833 C Check whether the specified bridging residues are cystines.
835 if (itype(iss(i)).ne.1) then
836 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
837 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
838 & ' can form a disulfide bridge?!!!'
839 write (*,'(2a,i3,a)')
840 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
841 & ' can form a disulfide bridge?!!!'
843 call MPI_Finalize(MPI_COMM_WORLD,ierror)
848 C Read preformed bridges.
850 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
851 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
854 C Check if the residues involved in bridges are in the specified list of
858 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
859 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
860 write (iout,'(a,i3,a)') 'Disulfide pair',i,
861 & ' contains residues present in other pairs.'
862 write (*,'(a,i3,a)') 'Disulfide pair',i,
863 & ' contains residues present in other pairs.'
865 call MPI_Finalize(MPI_COMM_WORLD,ierror)
871 if (ihpb(i).eq.iss(j)) goto 10
873 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
876 if (jhpb(i).eq.iss(j)) goto 20
878 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
891 c----------------------------------------------------------------------------
892 subroutine read_x(kanal,*)
893 implicit real*8 (a-h,o-z)
897 include 'COMMON.CHAIN'
898 include 'COMMON.IOUNITS'
899 include 'COMMON.CONTROL'
900 include 'COMMON.LOCAL'
901 include 'COMMON.INTERACT'
902 c Read coordinates from input
904 read(kanal,'(8f10.5)',end=10,err=10)
905 & ((c(l,k),l=1,3),k=1,nres),
906 & ((c(l,k+nres),l=1,3),k=nnt,nct)
909 c(j,2*nres)=c(j,nres)
911 call int_from_cart1(.false.)
914 dc(j,i)=c(j,i+1)-c(j,i)
915 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
919 if (itype(i).ne.10) then
921 dc(j,i+nres)=c(j,i+nres)-c(j,i)
922 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
930 c------------------------------------------------------------------------------
932 implicit real*8 (a-h,o-z)
934 include 'COMMON.IOUNITS'
937 include 'COMMON.INTERACT'
938 include 'COMMON.LOCAL'
939 include 'COMMON.NAMES'
940 include 'COMMON.CHAIN'
941 include 'COMMON.FFIELD'
942 include 'COMMON.SBRIDGE'
943 include 'COMMON.HEADER'
944 include 'COMMON.CONTROL'
945 c include 'COMMON.DBASE'
946 c include 'COMMON.THREAD'
947 include 'COMMON.TIME1'
948 C Set up variable list.
954 if (itype(i).ne.10) then
956 ialph(i,1)=nvar+nside
960 if (indphi.gt.0) then
962 else if (indback.gt.0) then
967 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
970 c----------------------------------------------------------------------------
971 subroutine gen_dist_constr
972 C Generate CA distance constraints.
973 implicit real*8 (a-h,o-z)
975 include 'COMMON.IOUNITS'
978 include 'COMMON.INTERACT'
979 include 'COMMON.LOCAL'
980 include 'COMMON.NAMES'
981 include 'COMMON.CHAIN'
982 include 'COMMON.FFIELD'
983 include 'COMMON.SBRIDGE'
984 include 'COMMON.HEADER'
985 include 'COMMON.CONTROL'
986 c include 'COMMON.DBASE'
987 c include 'COMMON.THREAD'
988 include 'COMMON.TIME1'
989 dimension itype_pdb(maxres)
990 common /pizda/ itype_pdb
992 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
993 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
994 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
996 do i=nstart_sup,nstart_sup+nsup-1
997 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
998 cd & ' seq_pdb', restyp(itype_pdb(i))
999 do j=i+2,nstart_sup+nsup-1
1001 ihpb(nhpb)=i+nstart_seq-nstart_sup
1002 jhpb(nhpb)=j+nstart_seq-nstart_sup
1004 dhpb(nhpb)=dist(i,j)
1007 cd write (iout,'(a)') 'Distance constraints:'
1012 cd if (ii.gt.nres) then
1017 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1018 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1019 cd & dhpb(i),forcon(i)
1024 subroutine read_minim
1025 implicit real*8 (a-h,o-z)
1026 include 'DIMENSIONS'
1027 include 'COMMON.MINIM'
1028 include 'COMMON.IOUNITS'
1030 character*320 minimcard
1031 call card_concat(minimcard)
1032 call readi(minimcard,'MAXMIN',maxmin,2000)
1033 call readi(minimcard,'MAXFUN',maxfun,5000)
1034 call readi(minimcard,'MINMIN',minmin,maxmin)
1035 call readi(minimcard,'MINFUN',minfun,maxmin)
1036 call reada(minimcard,'TOLF',tolf,1.0D-2)
1037 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1038 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1039 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1040 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1041 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1042 & 'Options in energy minimization:'
1043 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1044 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1045 & 'MinMin:',MinMin,' MinFun:',MinFun,
1046 & ' TolF:',TolF,' RTolF:',RTolF
1049 c----------------------------------------------------------------------------
1050 subroutine read_angles(kanal,*)
1051 implicit real*8 (a-h,o-z)
1052 include 'DIMENSIONS'
1053 include 'COMMON.GEO'
1054 include 'COMMON.VAR'
1055 include 'COMMON.CHAIN'
1056 include 'COMMON.IOUNITS'
1057 include 'COMMON.CONTROL'
1058 c Read angles from input
1060 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1061 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1062 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1063 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1066 c 9/7/01 avoid 180 deg valence angle
1067 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1069 theta(i)=deg2rad*theta(i)
1070 phi(i)=deg2rad*phi(i)
1071 alph(i)=deg2rad*alph(i)
1072 omeg(i)=deg2rad*omeg(i)
1077 c----------------------------------------------------------------------------
1078 subroutine reada(rekord,lancuch,wartosc,default)
1080 character*(*) rekord,lancuch
1081 double precision wartosc,default
1084 iread=index(rekord,lancuch)
1085 if (iread.eq.0) then
1089 iread=iread+ilen(lancuch)+1
1090 read (rekord(iread:),*,err=10,end=10) wartosc
1095 c----------------------------------------------------------------------------
1096 subroutine readi(rekord,lancuch,wartosc,default)
1098 character*(*) rekord,lancuch
1099 integer wartosc,default
1102 iread=index(rekord,lancuch)
1103 if (iread.eq.0) then
1107 iread=iread+ilen(lancuch)+1
1108 read (rekord(iread:),*,err=10,end=10) wartosc
1113 c----------------------------------------------------------------------------
1114 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1117 integer tablica(dim),default
1118 character*(*) rekord,lancuch
1125 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1126 if (iread.eq.0) return
1127 iread=iread+ilen(lancuch)+1
1128 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1131 c----------------------------------------------------------------------------
1132 subroutine multreada(rekord,lancuch,tablica,dim,default)
1135 double precision tablica(dim),default
1136 character*(*) rekord,lancuch
1143 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1144 if (iread.eq.0) return
1145 iread=iread+ilen(lancuch)+1
1146 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1149 c----------------------------------------------------------------------------
1150 subroutine openunits
1151 implicit real*8 (a-h,o-z)
1152 include 'DIMENSIONS'
1155 character*16 form,nodename
1158 include 'COMMON.SETUP'
1159 include 'COMMON.IOUNITS'
1160 c include 'COMMON.MD'
1161 include 'COMMON.CONTROL'
1162 integer lenpre,lenpot,ilen,lentmp
1164 character*3 out1file_text,ucase
1167 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1168 call getenv_loc("PREFIX",prefix)
1170 call getenv_loc("POT",pot)
1171 call getenv_loc("DIRTMP",tmpdir)
1172 call getenv_loc("CURDIR",curdir)
1173 call getenv_loc("OUT1FILE",out1file_text)
1174 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1175 out1file_text=ucase(out1file_text)
1176 if (out1file_text(1:1).eq."Y") then
1179 out1file=fg_rank.gt.0
1184 if (lentmp.gt.0) then
1185 write (*,'(80(1h!))')
1186 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1187 write (*,'(80(1h!))')
1188 write (*,*)"All output files will be on node /tmp directory."
1190 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1191 if (me.eq.king) then
1192 write (*,*) "The master node is ",nodename
1193 else if (fg_rank.eq.0) then
1194 write (*,*) "I am the CG slave node ",nodename
1196 write (*,*) "I am the FG slave node ",nodename
1199 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1200 lenpre = lentmp+lenpre+1
1202 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1203 C Get the names and open the input files
1204 #if defined(WINIFL) || defined(WINPGI)
1205 open(1,file=pref_orig(:ilen(pref_orig))//
1206 & '.inp',status='old',readonly,shared)
1207 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1208 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1209 C Get parameter filenames and open the parameter files.
1210 call getenv_loc('BONDPAR',bondname)
1211 open (ibond,file=bondname,status='old',readonly,shared)
1212 call getenv_loc('THETPAR',thetname)
1213 open (ithep,file=thetname,status='old',readonly,shared)
1215 call getenv_loc('THETPARPDB',thetname_pdb)
1216 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1218 call getenv_loc('ROTPAR',rotname)
1219 open (irotam,file=rotname,status='old',readonly,shared)
1221 call getenv_loc('ROTPARPDB',rotname_pdb)
1222 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1224 call getenv_loc('TORPAR',torname)
1225 open (itorp,file=torname,status='old',readonly,shared)
1226 call getenv_loc('TORDPAR',tordname)
1227 open (itordp,file=tordname,status='old',readonly,shared)
1228 call getenv_loc('FOURIER',fouriername)
1229 open (ifourier,file=fouriername,status='old',readonly,shared)
1230 call getenv_loc('ELEPAR',elename)
1231 open (ielep,file=elename,status='old',readonly,shared)
1232 call getenv_loc('SIDEPAR',sidename)
1233 open (isidep,file=sidename,status='old',readonly,shared)
1234 #elif (defined CRAY) || (defined AIX)
1235 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1237 c print *,"Processor",myrank," opened file 1"
1238 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1239 c print *,"Processor",myrank," opened file 9"
1240 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1241 C Get parameter filenames and open the parameter files.
1242 call getenv_loc('BONDPAR',bondname)
1243 open (ibond,file=bondname,status='old',action='read')
1244 c print *,"Processor",myrank," opened file IBOND"
1245 call getenv_loc('THETPAR',thetname)
1246 open (ithep,file=thetname,status='old',action='read')
1247 c print *,"Processor",myrank," opened file ITHEP"
1249 call getenv_loc('THETPARPDB',thetname_pdb)
1250 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1252 call getenv_loc('ROTPAR',rotname)
1253 open (irotam,file=rotname,status='old',action='read')
1254 c print *,"Processor",myrank," opened file IROTAM"
1256 call getenv_loc('ROTPARPDB',rotname_pdb)
1257 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1259 call getenv_loc('TORPAR',torname)
1260 open (itorp,file=torname,status='old',action='read')
1261 c print *,"Processor",myrank," opened file ITORP"
1262 call getenv_loc('TORDPAR',tordname)
1263 open (itordp,file=tordname,status='old',action='read')
1264 c print *,"Processor",myrank," opened file ITORDP"
1265 call getenv_loc('SCCORPAR',sccorname)
1266 open (isccor,file=sccorname,status='old',action='read')
1267 c print *,"Processor",myrank," opened file ISCCOR"
1268 call getenv_loc('FOURIER',fouriername)
1269 open (ifourier,file=fouriername,status='old',action='read')
1270 c print *,"Processor",myrank," opened file IFOURIER"
1271 call getenv_loc('ELEPAR',elename)
1272 open (ielep,file=elename,status='old',action='read')
1273 c print *,"Processor",myrank," opened file IELEP"
1274 call getenv_loc('SIDEPAR',sidename)
1275 open (isidep,file=sidename,status='old',action='read')
1276 c print *,"Processor",myrank," opened file ISIDEP"
1277 c print *,"Processor",myrank," opened parameter files"
1279 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1280 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1281 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1282 C Get parameter filenames and open the parameter files.
1283 call getenv_loc('BONDPAR',bondname)
1284 open (ibond,file=bondname,status='old')
1285 call getenv_loc('THETPAR',thetname)
1286 open (ithep,file=thetname,status='old')
1288 call getenv_loc('THETPARPDB',thetname_pdb)
1289 open (ithep_pdb,file=thetname_pdb,status='old')
1291 call getenv_loc('ROTPAR',rotname)
1292 open (irotam,file=rotname,status='old')
1294 call getenv_loc('ROTPARPDB',rotname_pdb)
1295 open (irotam_pdb,file=rotname_pdb,status='old')
1297 call getenv_loc('TORPAR',torname)
1298 open (itorp,file=torname,status='old')
1299 call getenv_loc('TORDPAR',tordname)
1300 open (itordp,file=tordname,status='old')
1301 call getenv_loc('SCCORPAR',sccorname)
1302 open (isccor,file=sccorname,status='old')
1303 call getenv_loc('FOURIER',fouriername)
1304 open (ifourier,file=fouriername,status='old')
1305 call getenv_loc('ELEPAR',elename)
1306 open (ielep,file=elename,status='old')
1307 call getenv_loc('SIDEPAR',sidename)
1308 open (isidep,file=sidename,status='old')
1310 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1312 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1313 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1314 C Get parameter filenames and open the parameter files.
1315 call getenv_loc('BONDPAR',bondname)
1316 open (ibond,file=bondname,status='old',readonly)
1317 call getenv_loc('THETPAR',thetname)
1318 open (ithep,file=thetname,status='old',readonly)
1320 call getenv_loc('THETPARPDB',thetname_pdb)
1321 print *,"thetname_pdb ",thetname_pdb
1322 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1323 print *,ithep_pdb," opened"
1325 call getenv_loc('ROTPAR',rotname)
1326 open (irotam,file=rotname,status='old',readonly)
1328 call getenv_loc('ROTPARPDB',rotname_pdb)
1329 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1331 call getenv_loc('TORPAR',torname)
1332 open (itorp,file=torname,status='old',readonly)
1333 call getenv_loc('TORDPAR',tordname)
1334 open (itordp,file=tordname,status='old',readonly)
1335 call getenv_loc('SCCORPAR',sccorname)
1336 open (isccor,file=sccorname,status='old',readonly)
1337 call getenv_loc('FOURIER',fouriername)
1338 open (ifourier,file=fouriername,status='old',readonly)
1339 call getenv_loc('ELEPAR',elename)
1340 open (ielep,file=elename,status='old',readonly)
1341 call getenv_loc('SIDEPAR',sidename)
1342 open (isidep,file=sidename,status='old',readonly)
1346 C 8/9/01 In the newest version SCp interaction constants are read from a file
1347 C Use -DOLDSCP to use hard-coded constants instead.
1349 call getenv_loc('SCPPAR',scpname)
1350 #if defined(WINIFL) || defined(WINPGI)
1351 open (iscpp,file=scpname,status='old',readonly,shared)
1352 #elif (defined CRAY) || (defined AIX)
1353 open (iscpp,file=scpname,status='old',action='read')
1355 open (iscpp,file=scpname,status='old')
1357 open (iscpp,file=scpname,status='old',readonly)
1360 call getenv_loc('PATTERN',patname)
1361 #if defined(WINIFL) || defined(WINPGI)
1362 open (icbase,file=patname,status='old',readonly,shared)
1363 #elif (defined CRAY) || (defined AIX)
1364 open (icbase,file=patname,status='old',action='read')
1366 open (icbase,file=patname,status='old')
1368 open (icbase,file=patname,status='old',readonly)
1371 C Open output file only for CG processes
1372 c print *,"Processor",myrank," fg_rank",fg_rank
1373 if (fg_rank.eq.0) then
1375 if (nodes.eq.1) then
1378 npos = dlog10(dfloat(nodes-1))+1
1380 if (npos.lt.3) npos=3
1381 write (liczba,'(i1)') npos
1382 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1384 write (liczba,form) me
1385 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1386 & liczba(:ilen(liczba))
1387 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1389 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1391 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1392 & liczba(:ilen(liczba))//'.mol2'
1393 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1394 & liczba(:ilen(liczba))//'.stat'
1396 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1397 & //liczba(:ilen(liczba))//'.stat')
1398 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1401 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1402 c & liczba(:ilen(liczba))//'.const'
1407 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1408 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1409 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1410 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1411 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1413 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1415 rest2name=prefix(:ilen(prefix))//'.rst'
1417 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1420 #if defined(AIX) || defined(PGI)
1421 if (me.eq.king .or. .not. out1file)
1422 & open(iout,file=outname,status='unknown')
1424 if (fg_rank.gt.0) then
1425 write (liczba,'(i3.3)') myrank/nfgtasks
1426 write (ll,'(bz,i3.3)') fg_rank
1427 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1432 open(igeom,file=intname,status='unknown',position='append')
1433 open(ipdb,file=pdbname,status='unknown')
1434 open(imol2,file=mol2name,status='unknown')
1435 open(istat,file=statname,status='unknown',position='append')
1437 c1out open(iout,file=outname,status='unknown')
1440 if (me.eq.king .or. .not.out1file)
1441 & open(iout,file=outname,status='unknown')
1443 if (fg_rank.gt.0) then
1444 write (liczba,'(i3.3)') myrank/nfgtasks
1445 write (ll,'(bz,i3.3)') fg_rank
1446 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1451 open(igeom,file=intname,status='unknown',access='append')
1452 open(ipdb,file=pdbname,status='unknown')
1453 open(imol2,file=mol2name,status='unknown')
1454 open(istat,file=statname,status='unknown',access='append')
1456 c1out open(iout,file=outname,status='unknown')
1459 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1460 csa_seed=prefix(:lenpre)//'.CSA.seed'
1461 csa_history=prefix(:lenpre)//'.CSA.history'
1462 csa_bank=prefix(:lenpre)//'.CSA.bank'
1463 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1464 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1465 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1466 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1467 csa_int=prefix(:lenpre)//'.int'
1468 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1469 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1470 csa_in=prefix(:lenpre)//'.CSA.in'
1471 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1474 write (iout,'(80(1h-))')
1475 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1476 write (iout,'(80(1h-))')
1477 write (iout,*) "Input file : ",
1478 & pref_orig(:ilen(pref_orig))//'.inp'
1479 write (iout,*) "Output file : ",
1480 & outname(:ilen(outname))
1482 write (iout,*) "Sidechain potential file : ",
1483 & sidename(:ilen(sidename))
1485 write (iout,*) "SCp potential file : ",
1486 & scpname(:ilen(scpname))
1488 write (iout,*) "Electrostatic potential file : ",
1489 & elename(:ilen(elename))
1490 write (iout,*) "Cumulant coefficient file : ",
1491 & fouriername(:ilen(fouriername))
1492 write (iout,*) "Torsional parameter file : ",
1493 & torname(:ilen(torname))
1494 write (iout,*) "Double torsional parameter file : ",
1495 & tordname(:ilen(tordname))
1496 write (iout,*) "SCCOR parameter file : ",
1497 & sccorname(:ilen(sccorname))
1498 write (iout,*) "Bond & inertia constant file : ",
1499 & bondname(:ilen(bondname))
1500 write (iout,*) "Bending parameter file : ",
1501 & thetname(:ilen(thetname))
1502 write (iout,*) "Rotamer parameter file : ",
1503 & rotname(:ilen(rotname))
1504 write (iout,*) "Threading database : ",
1505 & patname(:ilen(patname))
1507 &write (iout,*)" DIRTMP : ",
1509 write (iout,'(80(1h-))')
1513 c----------------------------------------------------------------------------
1514 subroutine card_concat(card)
1515 implicit real*8 (a-h,o-z)
1516 include 'DIMENSIONS'
1517 include 'COMMON.IOUNITS'
1519 character*80 karta,ucase
1521 read (inp,'(a)') karta
1524 do while (karta(80:80).eq.'&')
1525 card=card(:ilen(card)+1)//karta(:79)
1526 read (inp,'(a)') karta
1529 card=card(:ilen(card)+1)//karta
1533 subroutine read_dist_constr
1534 implicit real*8 (a-h,o-z)
1535 include 'DIMENSIONS'
1539 include 'COMMON.SETUP'
1540 include 'COMMON.CONTROL'
1541 include 'COMMON.CHAIN'
1542 include 'COMMON.IOUNITS'
1543 include 'COMMON.SBRIDGE'
1544 integer ifrag_(2,100),ipair_(2,100)
1545 double precision wfrag_(100),wpair_(100)
1546 character*500 controlcard
1547 c write (iout,*) "Calling read_dist_constr"
1548 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1550 call card_concat(controlcard)
1551 call readi(controlcard,"NFRAG",nfrag_,0)
1552 call readi(controlcard,"NPAIR",npair_,0)
1553 call readi(controlcard,"NDIST",ndist_,0)
1554 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1555 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1556 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1557 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1558 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1559 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1560 c write (iout,*) "IFRAG"
1562 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1564 c write (iout,*) "IPAIR"
1566 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1570 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1571 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1572 & ifrag_(2,i)=nstart_sup+nsup-1
1573 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1575 if (wfrag_(i).gt.0.0d0) then
1576 do j=ifrag_(1,i),ifrag_(2,i)-1
1577 do k=j+1,ifrag_(2,i)
1578 write (iout,*) "j",j," k",k
1580 if (constr_dist.eq.1) then
1585 forcon(nhpb)=wfrag_(i)
1586 else if (constr_dist.eq.2) then
1587 if (ddjk.le.dist_cut) then
1592 forcon(nhpb)=wfrag_(i)
1599 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1602 if (.not.out1file .or. me.eq.king)
1603 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1604 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1606 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1607 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1614 if (wpair_(i).gt.0.0d0) then
1622 do j=ifrag_(1,ii),ifrag_(2,ii)
1623 do k=ifrag_(1,jj),ifrag_(2,jj)
1627 forcon(nhpb)=wpair_(i)
1628 dhpb(nhpb)=dist(j,k)
1630 if (.not.out1file .or. me.eq.king)
1631 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1632 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1634 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1635 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1642 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1643 if (forcon(nhpb+1).gt.0.0d0) then
1645 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1647 if (.not.out1file .or. me.eq.king)
1648 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1649 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1651 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1652 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1659 c-------------------------------------------------------------------------------
1661 subroutine flush(iu)
1666 subroutine flush(iu)
1671 c------------------------------------------------------------------------------
1672 subroutine copy_to_tmp(source)
1673 include "DIMENSIONS"
1674 include "COMMON.IOUNITS"
1675 character*(*) source
1676 character* 256 tmpfile
1680 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1681 inquire(file=tmpfile,exist=ex)
1683 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1684 & " to temporary directory..."
1685 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1686 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1690 c------------------------------------------------------------------------------
1691 subroutine move_from_tmp(source)
1692 include "DIMENSIONS"
1693 include "COMMON.IOUNITS"
1694 character*(*) source
1697 write (*,*) "Moving ",source(:ilen(source)),
1698 & " from temporary directory to working directory"
1699 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1700 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1703 c------------------------------------------------------------------------------
1704 subroutine random_init(seed)
1706 C Initialize random number generator
1708 implicit real*8 (a-h,o-z)
1709 include 'DIMENSIONS'
1715 logical OKRandom, prng_restart
1717 integer iseed_array(4)
1719 include 'COMMON.IOUNITS'
1720 include 'COMMON.TIME1'
1721 c include 'COMMON.THREAD'
1722 include 'COMMON.SBRIDGE'
1723 include 'COMMON.CONTROL'
1724 include 'COMMON.MCM'
1725 c include 'COMMON.MAP'
1726 include 'COMMON.HEADER'
1727 c include 'COMMON.CSA'
1728 include 'COMMON.CHAIN'
1729 c include 'COMMON.MUCA'
1730 c include 'COMMON.MD'
1731 include 'COMMON.FFIELD'
1732 include 'COMMON.SETUP'
1733 iseed=-dint(dabs(seed))
1734 if (iseed.eq.0) then
1735 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1736 & 'Random seed undefined. The program will stop.'
1737 write (*,'(/80(1h*)/20x,a/80(1h*))')
1738 & 'Random seed undefined. The program will stop.'
1740 call mpi_finalize(mpi_comm_world,ierr)
1742 stop 'Bad random seed.'
1745 if (fg_rank.eq.0) then
1749 if(me.eq.king .or. .not. out1file)
1750 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1751 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1752 OKRandom = prng_restart(me,iseedi8)
1755 tmp=65536.0d0**(4-i)
1756 iseed_array(i) = dint(seed/tmp)
1757 seed=seed-iseed_array(i)*tmp
1759 if(me.eq.king .or. .not. out1file)
1760 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1761 & (iseed_array(i),i=1,4)
1762 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1763 & (iseed_array(i),i=1,4)
1764 OKRandom = prng_restart(me,iseed_array)
1767 c r1=ran_number(0.0D0,1.0D0)
1768 if(me.eq.king .or. .not. out1file)
1769 & write (iout,*) 'ran_num',r1
1770 if (r1.lt.0.0d0) OKRandom=.false.
1772 if (.not.OKRandom) then
1773 write (iout,*) 'PRNG IS NOT WORKING!!!'
1774 print *,'PRNG IS NOT WORKING!!!'
1777 call mpi_abort(mpi_comm_world,error_msg,ierr)
1780 write (iout,*) 'too many processors for parallel prng'
1781 write (*,*) 'too many processors for parallel prng'
1789 c write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)