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 if (modecalc.eq.8) then
32 inquire (file="fort.40",exist=file_exist)
33 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
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 c include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
74 c include 'COMMON.MAP'
75 include 'COMMON.HEADER'
76 c include 'COMMON.CSA'
77 include 'COMMON.CHAIN'
78 c include 'COMMON.MUCA'
80 include 'COMMON.FFIELD'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
100 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
101 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
102 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
103 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
104 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
105 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
106 call reada(controlcard,'DRMS',drms,0.1D0)
107 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
108 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
109 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
110 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
111 write (iout,'(a,f10.1)')'DRMS = ',drms
112 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
113 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
115 call readi(controlcard,'NZ_START',nz_start,0)
116 call readi(controlcard,'NZ_END',nz_end,0)
117 call readi(controlcard,'IZ_SC',iz_sc,0)
119 safety = 60.0d0*safety
122 call reada(controlcard,"T_BATH",t_bath,300.0d0)
123 minim=(index(controlcard,'MINIMIZE').gt.0)
124 dccart=(index(controlcard,'CART').gt.0)
125 overlapsc=(index(controlcard,'OVERLAP').gt.0)
126 overlapsc=.not.overlapsc
127 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
128 searchsc=.not.searchsc
129 sideadd=(index(controlcard,'SIDEADD').gt.0)
130 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
131 outpdb=(index(controlcard,'PDBOUT').gt.0)
132 outmol2=(index(controlcard,'MOL2OUT').gt.0)
133 pdbref=(index(controlcard,'PDBREF').gt.0)
134 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
135 indpdb=index(controlcard,'PDBSTART')
136 extconf=(index(controlcard,'EXTCONF').gt.0)
137 call readi(controlcard,'IPRINT',iprint,0)
138 call readi(controlcard,'MAXGEN',maxgen,10000)
139 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
140 call readi(controlcard,"KDIAG",kdiag,0)
141 call readi(controlcard,"RESCALE_MODE",rescale_mode,1)
142 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
143 & write (iout,*) "RESCALE_MODE",rescale_mode
144 split_ene=index(controlcard,'SPLIT_ENE').gt.0
145 if (index(controlcard,'REGULAR').gt.0.0D0) then
146 call reada(controlcard,'WEIDIS',weidis,0.1D0)
150 if (index(controlcard,'CHECKGRAD').gt.0) then
152 if (index(controlcard,'CART').gt.0) then
154 elseif (index(controlcard,'CARINT').gt.0) then
159 elseif (index(controlcard,'THREAD').gt.0) then
161 call readi(controlcard,'THREAD',nthread,0)
162 if (nthread.gt.0) then
163 call reada(controlcard,'WEIDIS',weidis,0.1D0)
166 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
167 stop 'Error termination in Read_Control.'
169 else if (index(controlcard,'MCMA').gt.0) then
171 else if (index(controlcard,'MCEE').gt.0) then
173 else if (index(controlcard,'MULTCONF').gt.0) then
175 else if (index(controlcard,'MAP').gt.0) then
177 call readi(controlcard,'MAP',nmap,0)
178 else if (index(controlcard,'CSA').gt.0) then
180 crc else if (index(controlcard,'ZSCORE').gt.0) then
182 crc ZSCORE is rm from UNRES, modecalc=9 is available
185 cfcm else if (index(controlcard,'MCMF').gt.0) then
187 else if (index(controlcard,'SOFTREG').gt.0) then
189 else if (index(controlcard,'CHECK_BOND').gt.0) then
191 else if (index(controlcard,'TEST').gt.0) then
193 else if (index(controlcard,'MD').gt.0) then
195 else if (index(controlcard,'RE ').gt.0) then
199 lmuca=index(controlcard,'MUCA').gt.0
200 call readi(controlcard,'MUCADYN',mucadyn,0)
201 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
202 if (lmuca .and. (me.eq.king .or. .not.out1file ))
204 write (iout,*) 'MUCADYN=',mucadyn
205 write (iout,*) 'MUCASMOOTH=',muca_smooth
208 iscode=index(controlcard,'ONE_LETTER')
209 indphi=index(controlcard,'PHI')
210 indback=index(controlcard,'BACK')
211 iranconf=index(controlcard,'RAND_CONF')
212 i2ndstr=index(controlcard,'USE_SEC_PRED')
213 gradout=index(controlcard,'GRADOUT').gt.0
214 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 if(me.eq.king.or..not.out1file)
217 & write (iout,'(2a)') diagmeth(kdiag),
218 & ' routine used to diagonalize matrices.'
221 c------------------------------------------------------------------------------
224 C Read molecular data.
226 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
235 include 'COMMON.INTERACT'
236 include 'COMMON.LOCAL'
237 include 'COMMON.NAMES'
238 include 'COMMON.CHAIN'
239 include 'COMMON.FFIELD'
240 include 'COMMON.SBRIDGE'
241 include 'COMMON.HEADER'
242 include 'COMMON.CONTROL'
243 c include 'COMMON.DBASE'
244 c include 'COMMON.THREAD'
245 include 'COMMON.CONTACTS'
246 include 'COMMON.TORCNSTR'
247 include 'COMMON.TIME1'
248 include 'COMMON.BOUNDS'
249 c include 'COMMON.MD'
250 c include 'COMMON.REMD'
251 include 'COMMON.SETUP'
252 character*4 sequence(maxres)
254 double precision x(maxvar)
255 character*256 pdbfile
256 character*320 weightcard
257 character*80 weightcard_t,ucase
258 dimension itype_pdb(maxres)
259 common /pizda/ itype_pdb
260 logical seq_comp,fail
261 double precision energia(0:n_ene)
267 C Read weights of the subsequent energy terms.
269 call card_concat(weightcard)
270 call reada(weightcard,'WLONG',wlong,1.0D0)
271 call reada(weightcard,'WSC',wsc,wlong)
272 call reada(weightcard,'WSCP',wscp,wlong)
273 call reada(weightcard,'WELEC',welec,1.0D0)
274 call reada(weightcard,'WVDWPP',wvdwpp,welec)
275 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
276 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
277 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
278 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
279 call reada(weightcard,'WTURN3',wturn3,1.0D0)
280 call reada(weightcard,'WTURN4',wturn4,1.0D0)
281 call reada(weightcard,'WTURN6',wturn6,1.0D0)
282 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
283 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
284 call reada(weightcard,'WBOND',wbond,1.0D0)
285 call reada(weightcard,'WTOR',wtor,1.0D0)
286 call reada(weightcard,'WTORD',wtor_d,1.0D0)
287 call reada(weightcard,'WANG',wang,1.0D0)
288 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
290 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
291 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
292 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
293 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
295 call reada(weightcard,'SCAL14',scal14,0.4D0)
296 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
297 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
298 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
299 call reada(weightcard,'TEMP0',temp0,300.0d0)
300 if (index(weightcard,'SOFT').gt.0) ipot=6
301 C 12/1/95 Added weight for the multi-body term WCORR
302 call reada(weightcard,'WCORRH',wcorr,1.0D0)
303 if (wcorr4.gt.0.0d0) wcorr=wcorr4
324 weights(24)=wdfa_dist
327 weights(27)=wdfa_beta
330 if(me.eq.king.or..not.out1file)
331 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
332 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
334 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
336 10 format (/'Energy-term weights (unscaled):'//
337 & 'WSCC= ',f10.6,' (SC-SC)'/
338 & 'WSCP= ',f10.6,' (SC-p)'/
339 & 'WELEC= ',f10.6,' (p-p electr)'/
340 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
341 & 'WBOND= ',f10.6,' (stretching)'/
342 & 'WANG= ',f10.6,' (bending)'/
343 & 'WSCLOC= ',f10.6,' (SC local)'/
344 & 'WTOR= ',f10.6,' (torsional)'/
345 & 'WTORD= ',f10.6,' (double torsional)'/
346 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
347 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
348 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
349 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
350 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
351 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
352 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
353 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
354 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
355 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
356 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
357 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
358 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
360 if(me.eq.king.or..not.out1file)then
361 if (wcorr4.gt.0.0d0) then
362 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
363 & 'between contact pairs of peptide groups'
364 write (iout,'(2(a,f5.3/))')
365 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
366 & 'Range of quenching the correlation terms:',2*delt_corr
367 else if (wcorr.gt.0.0d0) then
368 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
369 & 'between contact pairs of peptide groups'
371 write (iout,'(a,f8.3)')
372 & 'Scaling factor of 1,4 SC-p interactions:',scal14
373 write (iout,'(a,f8.3)')
374 & 'General scaling factor of SC-p interactions:',scalscp
376 r0_corr=cutoff_corr-delt_corr
378 aad(i,1)=scalscp*aad(i,1)
379 aad(i,2)=scalscp*aad(i,2)
380 bad(i,1)=scalscp*bad(i,1)
381 bad(i,2)=scalscp*bad(i,2)
383 call rescale_weights(t_bath)
384 if(me.eq.king.or..not.out1file)
385 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
386 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
388 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
390 22 format (/'Energy-term weights (scaled):'//
391 & 'WSCC= ',f10.6,' (SC-SC)'/
392 & 'WSCP= ',f10.6,' (SC-p)'/
393 & 'WELEC= ',f10.6,' (p-p electr)'/
394 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
395 & 'WBOND= ',f10.6,' (stretching)'/
396 & 'WANG= ',f10.6,' (bending)'/
397 & 'WSCLOC= ',f10.6,' (SC local)'/
398 & 'WTOR= ',f10.6,' (torsional)'/
399 & 'WTORD= ',f10.6,' (double torsional)'/
400 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
401 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
402 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
403 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
404 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
405 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
406 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
407 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
408 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
409 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
410 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
411 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
412 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
414 if(me.eq.king.or..not.out1file)
415 & write (iout,*) "Reference temperature for weights calculation:",
417 call reada(weightcard,"D0CM",d0cm,3.78d0)
418 call reada(weightcard,"AKCM",akcm,15.1d0)
419 call reada(weightcard,"AKTH",akth,11.0d0)
420 call reada(weightcard,"AKCT",akct,12.0d0)
421 call reada(weightcard,"V1SS",v1ss,-1.08d0)
422 call reada(weightcard,"V2SS",v2ss,7.61d0)
423 call reada(weightcard,"V3SS",v3ss,13.7d0)
424 call reada(weightcard,"EBR",ebr,-5.50D0)
425 if(me.eq.king.or..not.out1file) then
426 write (iout,*) "Parameters of the SS-bond potential:"
427 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
429 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
430 write (iout,*) "EBR",ebr
431 print *,'indpdb=',indpdb,' pdbref=',pdbref
433 if (indpdb.gt.0 .or. pdbref) then
434 read(inp,'(a)') pdbfile
435 if(me.eq.king.or..not.out1file)
436 & write (iout,'(2a)') 'PDB data will be read from file ',
437 & pdbfile(:ilen(pdbfile))
438 open(ipdbin,file=pdbfile,status='old',err=33)
440 33 write (iout,'(a)') 'Error opening PDB file.'
443 c print *,'Begin reading pdb data'
445 c print *,'Finished reading pdb data'
446 if(me.eq.king.or..not.out1file)
447 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
448 & ' nstart_sup=',nstart_sup
450 itype_pdb(i)=itype(i)
454 nct=nstart_sup+nsup-1
455 call contact(.false.,ncont_ref,icont_ref,co)
458 if(me.eq.king.or..not.out1file)
459 & write(iout,*)'Adding sidechains'
466 do while (fail.and.nsi.le.maxsi)
467 c call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
470 if(fail) write(iout,*)'Adding sidechain failed for res ',
471 & i,' after ',nsi,' trials'
477 if (indpdb.eq.0) then
478 C Read sequence if not taken from the pdb file.
480 c print *,'nres=',nres
481 if (iscode.gt.0) then
482 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
484 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
486 C Convert sequence to numeric code
488 itype(i)=rescode(i,sequence(i),iscode)
490 C Assign initial virtual bond lengths
496 vbld(i+nres)=dsc(itype(i))
497 vbld_inv(i+nres)=dsc_inv(itype(i))
498 c write (iout,*) "i",i," itype",itype(i),
499 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
503 c print '(20i4)',(itype(i),i=1,nres)
506 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
508 if (itype(i).eq.21) then
512 else if (itype(i+1).ne.20) then
514 else if (itype(i).ne.20) then
521 if(me.eq.king.or..not.out1file)then
522 write (iout,*) "ITEL"
524 write (iout,*) i,itype(i),itel(i)
526 print *,'Call Read_Bridge.'
529 C 8/13/98 Set limits to generating the dihedral angles
534 read (inp,*) ndih_constr
535 if (ndih_constr.gt.0) then
537 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
538 if(me.eq.king.or..not.out1file)then
540 & 'There are',ndih_constr,' constraints on phi angles.'
542 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
546 phi0(i)=deg2rad*phi0(i)
547 drange(i)=deg2rad*drange(i)
549 if(me.eq.king.or..not.out1file)
550 & write (iout,*) 'FTORS',ftors
553 phibound(1,ii) = phi0(i)-drange(i)
554 phibound(2,ii) = phi0(i)+drange(i)
561 write (iout,'(a)') 'Boundaries in phi angle sampling:'
563 write (iout,'(a3,i5,2f10.1)')
564 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
570 cd print *,'NNT=',NNT,' NCT=',NCT
571 if (itype(1).eq.21) nnt=2
572 if (itype(nres).eq.21) nct=nct-1
574 C Juyong:READ init_vars
575 C Initialize variables!
576 C Juyong:READ read_info
577 C READ fragment information!!
578 C both routines should be in dfa.F file!!
580 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
581 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
583 print*, 'init_dfa_vars finished!'
585 print*, 'read_dfa_info finished!'
592 if(me.eq.king.or..not.out1file)
593 & write (iout,'(a,i3)') 'nsup=',nsup
595 if (nsup.le.(nct-nnt+1)) then
596 do i=0,nct-nnt+1-nsup
597 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
603 & 'Error - sequences to be superposed do not match.'
606 do i=0,nsup-(nct-nnt+1)
607 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
609 nstart_sup=nstart_sup+i
615 & 'Error - sequences to be superposed do not match.'
618 if (nsup.eq.0) nsup=nct-nnt
619 if (nstart_sup.eq.0) nstart_sup=nnt
620 if (nstart_seq.eq.0) nstart_seq=nnt
621 if(me.eq.king.or..not.out1file)
622 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
623 & ' nstart_seq=',nstart_seq
625 c--- Zscore rms -------
626 if (nz_start.eq.0) nz_start=nnt
627 if (nz_end.eq.0 .and. nsup.gt.0) then
629 else if (nz_end.eq.0) then
632 if(me.eq.king.or..not.out1file)then
633 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
634 write (iout,*) 'IZ_SC=',iz_sc
636 c----------------------
639 if (.not.pdbref) then
640 call read_angles(inp,*38)
642 38 write (iout,'(a)') 'Error reading reference structure.'
644 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
645 stop 'Error reading reference structure'
649 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
658 call contact(.true.,ncont_ref,icont_ref,co)
660 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
662 if (constr_dist.gt.0) call read_dist_constr
663 c write (iout,*) "After read_dist_constr nhpb",nhpb
665 if(me.eq.king.or..not.out1file)
666 & write (iout,*) 'Contact order:',co
668 if(me.eq.king.or..not.out1file)
669 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
672 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
674 if(me.eq.king.or..not.out1file)
675 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
676 & icont_ref(1,i),' ',
677 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
681 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
682 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
683 & modecalc.ne.10) then
684 C If input structure hasn't been supplied from the PDB file read or generate
686 if (iranconf.eq.0 .and. .not. extconf) then
687 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
688 & write (iout,'(a)') 'Initial geometry will be read in.'
690 read(inp,'(8f10.5)',end=36,err=36)
691 & ((c(l,k),l=1,3),k=1,nres),
692 & ((c(l,k+nres),l=1,3),k=nnt,nct)
693 call int_from_cart1(.false.)
696 dc(j,i)=c(j,i+1)-c(j,i)
697 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
701 if (itype(i).ne.10) then
703 dc(j,i+nres)=c(j,i+nres)-c(j,i)
704 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
710 call read_angles(inp,*36)
713 36 write (iout,'(a)') 'Error reading angle file.'
715 call mpi_finalize( MPI_COMM_WORLD,IERR )
717 stop 'Error reading angle file.'
719 else if (extconf) then
720 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
721 & write (iout,'(a)') 'Extended chain initial geometry.'
723 theta(i)=90d0*deg2rad
729 alph(i)=110d0*deg2rad
732 omeg(i)=-120d0*deg2rad
735 if(me.eq.king.or..not.out1file)
736 & write (iout,'(a)') 'Random-generated initial geometry.'
740 if (me.eq.king .or. fg_rank.eq.0 .and. (
741 & modecalc.eq.12 .or. modecalc.eq.14) ) then
745 c call gen_rand_conf(itmp,*30)
747 30 write (iout,*) 'Failed to generate random conformation',
749 write (*,*) 'Processor:',me,
750 & ' Failed to generate random conformation',
760 write (iout,'(a,i3,a)') 'Processor:',me,
761 & ' error in generating random conformation.'
762 write (*,'(a,i3,a)') 'Processor:',me,
763 & ' error in generating random conformation.'
766 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
772 c call gen_rand_conf(itmp,*31)
774 31 write (iout,*) 'Failed to generate random conformation',
776 write (*,*) 'Failed to generate random conformation',
779 write (iout,'(a,i3,a)') 'Processor:',me,
780 & ' error in generating random conformation.'
781 write (*,'(a,i3,a)') 'Processor:',me,
782 & ' error in generating random conformation.'
787 elseif (modecalc.eq.4) then
788 read (inp,'(a)') intinname
789 open (intin,file=intinname,status='old',err=333)
790 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
791 & write (iout,'(a)') 'intinname',intinname
792 write (*,'(a)') 'Processor',myrank,' intinname',intinname
794 333 write (iout,'(2a)') 'Error opening angle file ',intinname
796 call MPI_Finalize(MPI_COMM_WORLD,IERR)
798 stop 'Error opening angle file.'
802 C Generate distance constraints, if the PDB structure is to be regularized.
803 if (nthread.gt.0) then
807 if (me.eq.king .or. .not. out1file)
809 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
810 write (iout,'(/a,i3,a)')
811 & 'The chain contains',ns,' disulfide-bridging cysteines.'
812 write (iout,'(20i4)') (iss(i),i=1,ns)
813 write (iout,'(/a/)') 'Pre-formed links are:'
819 if (me.eq.king.or..not.out1file)
820 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
821 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
826 c if (i2ndstr.gt.0) call secstrp2dihc
827 c call geom_to_var(nvar,x)
828 c call etotal(energia(0))
829 c call enerprint(energia(0))
830 c call briefout(0,etot)
832 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
833 cd write (iout,'(a)') 'Variable list:'
834 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
836 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
837 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
838 & 'Processor',myrank,': end reading molecular data.'
842 c--------------------------------------------------------------------------
843 logical function seq_comp(itypea,itypeb,length)
845 integer length,itypea(length),itypeb(length)
848 if (itypea(i).ne.itypeb(i)) then
856 c-----------------------------------------------------------------------------
857 subroutine read_bridge
858 C Read information about disulfide bridges.
859 implicit real*8 (a-h,o-z)
864 include 'COMMON.IOUNITS'
867 include 'COMMON.INTERACT'
868 include 'COMMON.LOCAL'
869 include 'COMMON.NAMES'
870 include 'COMMON.CHAIN'
871 include 'COMMON.FFIELD'
872 include 'COMMON.SBRIDGE'
873 include 'COMMON.HEADER'
874 include 'COMMON.CONTROL'
875 c include 'COMMON.DBASE'
876 c include 'COMMON.THREAD'
877 include 'COMMON.TIME1'
878 include 'COMMON.SETUP'
879 C Read bridging residues.
880 read (inp,*) ns,(iss(i),i=1,ns)
882 if(me.eq.king.or..not.out1file)
883 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
884 C Check whether the specified bridging residues are cystines.
886 if (itype(iss(i)).ne.1) then
887 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
888 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
889 & ' can form a disulfide bridge?!!!'
890 write (*,'(2a,i3,a)')
891 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
892 & ' can form a disulfide bridge?!!!'
894 call MPI_Finalize(MPI_COMM_WORLD,ierror)
899 C Read preformed bridges.
901 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
902 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
905 C Check if the residues involved in bridges are in the specified list of
909 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
910 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
911 write (iout,'(a,i3,a)') 'Disulfide pair',i,
912 & ' contains residues present in other pairs.'
913 write (*,'(a,i3,a)') 'Disulfide pair',i,
914 & ' contains residues present in other pairs.'
916 call MPI_Finalize(MPI_COMM_WORLD,ierror)
922 if (ihpb(i).eq.iss(j)) goto 10
924 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
927 if (jhpb(i).eq.iss(j)) goto 20
929 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
942 c----------------------------------------------------------------------------
943 subroutine read_x(kanal,*)
944 implicit real*8 (a-h,o-z)
948 include 'COMMON.CHAIN'
949 include 'COMMON.IOUNITS'
950 include 'COMMON.CONTROL'
951 include 'COMMON.LOCAL'
952 include 'COMMON.INTERACT'
953 c Read coordinates from input
955 read(kanal,'(8f10.5)',end=10,err=10)
956 & ((c(l,k),l=1,3),k=1,nres),
957 & ((c(l,k+nres),l=1,3),k=nnt,nct)
960 c(j,2*nres)=c(j,nres)
962 call int_from_cart1(.false.)
965 dc(j,i)=c(j,i+1)-c(j,i)
966 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
970 if (itype(i).ne.10) then
972 dc(j,i+nres)=c(j,i+nres)-c(j,i)
973 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
981 c------------------------------------------------------------------------------
983 implicit real*8 (a-h,o-z)
985 include 'COMMON.IOUNITS'
988 include 'COMMON.INTERACT'
989 include 'COMMON.LOCAL'
990 include 'COMMON.NAMES'
991 include 'COMMON.CHAIN'
992 include 'COMMON.FFIELD'
993 include 'COMMON.SBRIDGE'
994 include 'COMMON.HEADER'
995 include 'COMMON.CONTROL'
996 c include 'COMMON.DBASE'
997 c include 'COMMON.THREAD'
998 include 'COMMON.TIME1'
999 C Set up variable list.
1005 if (itype(i).ne.10) then
1007 ialph(i,1)=nvar+nside
1011 if (indphi.gt.0) then
1013 else if (indback.gt.0) then
1018 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1021 c----------------------------------------------------------------------------
1022 subroutine gen_dist_constr
1023 C Generate CA distance constraints.
1024 implicit real*8 (a-h,o-z)
1025 include 'DIMENSIONS'
1026 include 'COMMON.IOUNITS'
1027 include 'COMMON.GEO'
1028 include 'COMMON.VAR'
1029 include 'COMMON.INTERACT'
1030 include 'COMMON.LOCAL'
1031 include 'COMMON.NAMES'
1032 include 'COMMON.CHAIN'
1033 include 'COMMON.FFIELD'
1034 include 'COMMON.SBRIDGE'
1035 include 'COMMON.HEADER'
1036 include 'COMMON.CONTROL'
1037 c include 'COMMON.DBASE'
1038 c include 'COMMON.THREAD'
1039 include 'COMMON.TIME1'
1040 dimension itype_pdb(maxres)
1041 common /pizda/ itype_pdb
1043 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1044 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1045 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1047 do i=nstart_sup,nstart_sup+nsup-1
1048 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1049 cd & ' seq_pdb', restyp(itype_pdb(i))
1050 do j=i+2,nstart_sup+nsup-1
1052 ihpb(nhpb)=i+nstart_seq-nstart_sup
1053 jhpb(nhpb)=j+nstart_seq-nstart_sup
1055 dhpb(nhpb)=dist(i,j)
1058 cd write (iout,'(a)') 'Distance constraints:'
1063 cd if (ii.gt.nres) then
1068 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1069 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1070 cd & dhpb(i),forcon(i)
1074 c----------------------------------------------------------------------------
1076 implicit real*8 (a-h,o-z)
1077 include 'DIMENSIONS'
1078 include 'COMMON.IOUNITS'
1079 include 'COMMON.GEO'
1080 include 'COMMON.CSA'
1081 include 'COMMON.BANK'
1082 include 'COMMON.CONTROL'
1084 character*620 mcmcard
1085 call card_concat(mcmcard)
1087 call readi(mcmcard,'NCONF',nconf,50)
1088 call readi(mcmcard,'NADD',nadd,0)
1089 call readi(mcmcard,'JSTART',jstart,1)
1090 call readi(mcmcard,'JEND',jend,1)
1091 call readi(mcmcard,'NSTMAX',nstmax,500000)
1092 call readi(mcmcard,'N0',n0,1)
1093 call readi(mcmcard,'N1',n1,6)
1094 call readi(mcmcard,'N2',n2,4)
1095 call readi(mcmcard,'N3',n3,0)
1096 call readi(mcmcard,'N4',n4,0)
1097 call readi(mcmcard,'N5',n5,0)
1098 call readi(mcmcard,'N6',n6,10)
1099 call readi(mcmcard,'N7',n7,0)
1100 call readi(mcmcard,'N8',n8,0)
1101 call readi(mcmcard,'N9',n9,0)
1102 call readi(mcmcard,'N14',n14,0)
1103 call readi(mcmcard,'N15',n15,0)
1104 call readi(mcmcard,'N16',n16,0)
1105 call readi(mcmcard,'N17',n17,0)
1106 call readi(mcmcard,'N18',n18,0)
1108 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1110 call readi(mcmcard,'NDIFF',ndiff,2)
1111 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1112 call readi(mcmcard,'IS1',is1,1)
1113 call readi(mcmcard,'IS2',is2,8)
1114 call readi(mcmcard,'NRAN0',nran0,4)
1115 call readi(mcmcard,'NRAN1',nran1,2)
1116 call readi(mcmcard,'IRR',irr,1)
1117 call readi(mcmcard,'NSEED',nseed,20)
1118 call readi(mcmcard,'NTOTAL',ntotal,10000)
1119 call reada(mcmcard,'CUT1',cut1,2.0d0)
1120 call reada(mcmcard,'CUT2',cut2,5.0d0)
1121 call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1122 call readi(mcmcard,'ICMAX',icmax,1)
1123 call readi(mcmcard,'NBANKM',nbankm,400)
1124 call readi(mcmcard,'IUCUT',iucut,2)
1125 call readi(mcmcard,'IRESTART',irestart,0)
1126 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1129 call reada(mcmcard,'DELE',dele,20.0d0)
1130 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1131 call readi(mcmcard,'IREF',iref,0)
1132 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1133 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1134 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1135 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1136 write (iout,*) "NCONF_IN",nconf_in
1137 tm_score=(index(mcmcard,'TMSCORE').gt.0)
1138 if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1139 & " for torsional angles"
1143 subroutine read_minim
1144 implicit real*8 (a-h,o-z)
1145 include 'DIMENSIONS'
1146 include 'COMMON.MINIM'
1147 include 'COMMON.IOUNITS'
1149 character*320 minimcard
1150 call card_concat(minimcard)
1151 call readi(minimcard,'MAXMIN',maxmin,2000)
1152 call readi(minimcard,'MAXFUN',maxfun,5000)
1153 call readi(minimcard,'MINMIN',minmin,maxmin)
1154 call readi(minimcard,'MINFUN',minfun,maxmin)
1155 call reada(minimcard,'TOLF',tolf,1.0D-2)
1156 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1157 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1158 & 'Options in energy minimization:'
1159 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1160 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1161 & 'MinMin:',MinMin,' MinFun:',MinFun,
1162 & ' TolF:',TolF,' RTolF:',RTolF
1165 c----------------------------------------------------------------------------
1166 subroutine read_angles(kanal,*)
1167 implicit real*8 (a-h,o-z)
1168 include 'DIMENSIONS'
1169 include 'COMMON.GEO'
1170 include 'COMMON.VAR'
1171 include 'COMMON.CHAIN'
1172 include 'COMMON.IOUNITS'
1173 include 'COMMON.CONTROL'
1174 c Read angles from input
1176 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1177 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1178 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1179 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1182 c 9/7/01 avoid 180 deg valence angle
1183 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1185 theta(i)=deg2rad*theta(i)
1186 phi(i)=deg2rad*phi(i)
1187 alph(i)=deg2rad*alph(i)
1188 omeg(i)=deg2rad*omeg(i)
1193 c----------------------------------------------------------------------------
1194 subroutine reada(rekord,lancuch,wartosc,default)
1196 character*(*) rekord,lancuch
1197 double precision wartosc,default
1200 iread=index(rekord,lancuch)
1201 if (iread.eq.0) then
1205 iread=iread+ilen(lancuch)+1
1206 read (rekord(iread:),*,err=10,end=10) wartosc
1211 c----------------------------------------------------------------------------
1212 subroutine readi(rekord,lancuch,wartosc,default)
1214 character*(*) rekord,lancuch
1215 integer wartosc,default
1218 iread=index(rekord,lancuch)
1219 if (iread.eq.0) then
1223 iread=iread+ilen(lancuch)+1
1224 read (rekord(iread:),*,err=10,end=10) wartosc
1229 c----------------------------------------------------------------------------
1230 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1233 integer tablica(dim),default
1234 character*(*) rekord,lancuch
1241 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1242 if (iread.eq.0) return
1243 iread=iread+ilen(lancuch)+1
1244 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1247 c----------------------------------------------------------------------------
1248 subroutine multreada(rekord,lancuch,tablica,dim,default)
1251 double precision tablica(dim),default
1252 character*(*) rekord,lancuch
1259 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1260 if (iread.eq.0) return
1261 iread=iread+ilen(lancuch)+1
1262 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1265 c----------------------------------------------------------------------------
1266 subroutine openunits
1267 implicit real*8 (a-h,o-z)
1268 include 'DIMENSIONS'
1271 character*16 form,nodename
1274 include 'COMMON.SETUP'
1275 include 'COMMON.IOUNITS'
1276 c include 'COMMON.MD'
1277 include 'COMMON.CONTROL'
1278 integer lenpre,lenpot,ilen,lentmp
1280 character*3 out1file_text,ucase
1283 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1284 call getenv_loc("PREFIX",prefix)
1286 call getenv_loc("POT",pot)
1287 call getenv_loc("DIRTMP",tmpdir)
1288 call getenv_loc("CURDIR",curdir)
1289 call getenv_loc("OUT1FILE",out1file_text)
1290 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1291 out1file_text=ucase(out1file_text)
1292 if (out1file_text(1:1).eq."Y") then
1295 out1file=fg_rank.gt.0
1300 if (lentmp.gt.0) then
1301 write (*,'(80(1h!))')
1302 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1303 write (*,'(80(1h!))')
1304 write (*,*)"All output files will be on node /tmp directory."
1306 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1307 if (me.eq.king) then
1308 write (*,*) "The master node is ",nodename
1309 else if (fg_rank.eq.0) then
1310 write (*,*) "I am the CG slave node ",nodename
1312 write (*,*) "I am the FG slave node ",nodename
1315 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1316 lenpre = lentmp+lenpre+1
1318 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1319 C Get the names and open the input files
1320 #if defined(WINIFL) || defined(WINPGI)
1321 open(1,file=pref_orig(:ilen(pref_orig))//
1322 & '.inp',status='old',readonly,shared)
1323 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1324 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1325 C Get parameter filenames and open the parameter files.
1326 call getenv_loc('BONDPAR',bondname)
1327 open (ibond,file=bondname,status='old',readonly,shared)
1328 call getenv_loc('THETPAR',thetname)
1329 open (ithep,file=thetname,status='old',readonly,shared)
1331 call getenv_loc('THETPARPDB',thetname_pdb)
1332 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1334 call getenv_loc('ROTPAR',rotname)
1335 open (irotam,file=rotname,status='old',readonly,shared)
1337 call getenv_loc('ROTPARPDB',rotname_pdb)
1338 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1340 call getenv_loc('TORPAR',torname)
1341 open (itorp,file=torname,status='old',readonly,shared)
1342 call getenv_loc('TORDPAR',tordname)
1343 open (itordp,file=tordname,status='old',readonly,shared)
1344 call getenv_loc('FOURIER',fouriername)
1345 open (ifourier,file=fouriername,status='old',readonly,shared)
1346 call getenv_loc('ELEPAR',elename)
1347 open (ielep,file=elename,status='old',readonly,shared)
1348 call getenv_loc('SIDEPAR',sidename)
1349 open (isidep,file=sidename,status='old',readonly,shared)
1350 #elif (defined CRAY) || (defined AIX)
1351 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1353 c print *,"Processor",myrank," opened file 1"
1354 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1355 c print *,"Processor",myrank," opened file 9"
1356 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1357 C Get parameter filenames and open the parameter files.
1358 call getenv_loc('BONDPAR',bondname)
1359 open (ibond,file=bondname,status='old',action='read')
1360 c print *,"Processor",myrank," opened file IBOND"
1361 call getenv_loc('THETPAR',thetname)
1362 open (ithep,file=thetname,status='old',action='read')
1363 c print *,"Processor",myrank," opened file ITHEP"
1365 call getenv_loc('THETPARPDB',thetname_pdb)
1366 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1368 call getenv_loc('ROTPAR',rotname)
1369 open (irotam,file=rotname,status='old',action='read')
1370 c print *,"Processor",myrank," opened file IROTAM"
1372 call getenv_loc('ROTPARPDB',rotname_pdb)
1373 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1375 call getenv_loc('TORPAR',torname)
1376 open (itorp,file=torname,status='old',action='read')
1377 c print *,"Processor",myrank," opened file ITORP"
1378 call getenv_loc('TORDPAR',tordname)
1379 open (itordp,file=tordname,status='old',action='read')
1380 c print *,"Processor",myrank," opened file ITORDP"
1381 call getenv_loc('SCCORPAR',sccorname)
1382 open (isccor,file=sccorname,status='old',action='read')
1383 c print *,"Processor",myrank," opened file ISCCOR"
1384 call getenv_loc('FOURIER',fouriername)
1385 open (ifourier,file=fouriername,status='old',action='read')
1386 c print *,"Processor",myrank," opened file IFOURIER"
1387 call getenv_loc('ELEPAR',elename)
1388 open (ielep,file=elename,status='old',action='read')
1389 c print *,"Processor",myrank," opened file IELEP"
1390 call getenv_loc('SIDEPAR',sidename)
1391 open (isidep,file=sidename,status='old',action='read')
1392 c print *,"Processor",myrank," opened file ISIDEP"
1393 c print *,"Processor",myrank," opened parameter files"
1395 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1396 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1397 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1398 C Get parameter filenames and open the parameter files.
1399 call getenv_loc('BONDPAR',bondname)
1400 open (ibond,file=bondname,status='old')
1401 call getenv_loc('THETPAR',thetname)
1402 open (ithep,file=thetname,status='old')
1404 call getenv_loc('THETPARPDB',thetname_pdb)
1405 open (ithep_pdb,file=thetname_pdb,status='old')
1407 call getenv_loc('ROTPAR',rotname)
1408 open (irotam,file=rotname,status='old')
1410 call getenv_loc('ROTPARPDB',rotname_pdb)
1411 open (irotam_pdb,file=rotname_pdb,status='old')
1413 call getenv_loc('TORPAR',torname)
1414 open (itorp,file=torname,status='old')
1415 call getenv_loc('TORDPAR',tordname)
1416 open (itordp,file=tordname,status='old')
1417 call getenv_loc('SCCORPAR',sccorname)
1418 open (isccor,file=sccorname,status='old')
1419 call getenv_loc('FOURIER',fouriername)
1420 open (ifourier,file=fouriername,status='old')
1421 call getenv_loc('ELEPAR',elename)
1422 open (ielep,file=elename,status='old')
1423 call getenv_loc('SIDEPAR',sidename)
1424 open (isidep,file=sidename,status='old')
1426 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1428 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1429 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1430 C Get parameter filenames and open the parameter files.
1431 call getenv_loc('BONDPAR',bondname)
1432 open (ibond,file=bondname,status='old',readonly)
1433 call getenv_loc('THETPAR',thetname)
1434 open (ithep,file=thetname,status='old',readonly)
1436 call getenv_loc('THETPARPDB',thetname_pdb)
1437 print *,"thetname_pdb ",thetname_pdb
1438 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1439 print *,ithep_pdb," opened"
1441 call getenv_loc('ROTPAR',rotname)
1442 open (irotam,file=rotname,status='old',readonly)
1444 call getenv_loc('ROTPARPDB',rotname_pdb)
1445 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1447 call getenv_loc('TORPAR',torname)
1448 open (itorp,file=torname,status='old',readonly)
1449 call getenv_loc('TORDPAR',tordname)
1450 open (itordp,file=tordname,status='old',readonly)
1451 call getenv_loc('SCCORPAR',sccorname)
1452 open (isccor,file=sccorname,status='old',readonly)
1453 call getenv_loc('FOURIER',fouriername)
1454 open (ifourier,file=fouriername,status='old',readonly)
1455 call getenv_loc('ELEPAR',elename)
1456 open (ielep,file=elename,status='old',readonly)
1457 call getenv_loc('SIDEPAR',sidename)
1458 open (isidep,file=sidename,status='old',readonly)
1462 C 8/9/01 In the newest version SCp interaction constants are read from a file
1463 C Use -DOLDSCP to use hard-coded constants instead.
1465 call getenv_loc('SCPPAR',scpname)
1466 #if defined(WINIFL) || defined(WINPGI)
1467 open (iscpp,file=scpname,status='old',readonly,shared)
1468 #elif (defined CRAY) || (defined AIX)
1469 open (iscpp,file=scpname,status='old',action='read')
1471 open (iscpp,file=scpname,status='old')
1473 open (iscpp,file=scpname,status='old',readonly)
1476 call getenv_loc('PATTERN',patname)
1477 #if defined(WINIFL) || defined(WINPGI)
1478 open (icbase,file=patname,status='old',readonly,shared)
1479 #elif (defined CRAY) || (defined AIX)
1480 open (icbase,file=patname,status='old',action='read')
1482 open (icbase,file=patname,status='old')
1484 open (icbase,file=patname,status='old',readonly)
1487 C Open output file only for CG processes
1488 c print *,"Processor",myrank," fg_rank",fg_rank
1489 if (fg_rank.eq.0) then
1491 if (nodes.eq.1) then
1494 npos = dlog10(dfloat(nodes-1))+1
1496 if (npos.lt.3) npos=3
1497 write (liczba,'(i1)') npos
1498 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1500 write (liczba,form) me
1501 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1502 & liczba(:ilen(liczba))
1503 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1505 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1507 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1508 & liczba(:ilen(liczba))//'.mol2'
1509 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1510 & liczba(:ilen(liczba))//'.stat'
1512 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1513 & //liczba(:ilen(liczba))//'.stat')
1514 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1517 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1518 c & liczba(:ilen(liczba))//'.const'
1523 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1524 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1525 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1526 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1527 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1529 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1531 rest2name=prefix(:ilen(prefix))//'.rst'
1533 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1536 #if defined(AIX) || defined(PGI)
1537 if (me.eq.king .or. .not. out1file)
1538 & open(iout,file=outname,status='unknown')
1540 if (fg_rank.gt.0) then
1541 write (liczba,'(i3.3)') myrank/nfgtasks
1542 write (ll,'(bz,i3.3)') fg_rank
1543 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1548 open(igeom,file=intname,status='unknown',position='append')
1549 open(ipdb,file=pdbname,status='unknown')
1550 open(imol2,file=mol2name,status='unknown')
1551 open(istat,file=statname,status='unknown',position='append')
1553 c1out open(iout,file=outname,status='unknown')
1556 if (me.eq.king .or. .not.out1file)
1557 & open(iout,file=outname,status='unknown')
1559 if (fg_rank.gt.0) then
1560 write (liczba,'(i3.3)') myrank/nfgtasks
1561 write (ll,'(bz,i3.3)') fg_rank
1562 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1567 open(igeom,file=intname,status='unknown',access='append')
1568 open(ipdb,file=pdbname,status='unknown')
1569 open(imol2,file=mol2name,status='unknown')
1570 open(istat,file=statname,status='unknown',access='append')
1572 c1out open(iout,file=outname,status='unknown')
1575 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1576 csa_seed=prefix(:lenpre)//'.CSA.seed'
1577 csa_history=prefix(:lenpre)//'.CSA.history'
1578 csa_bank=prefix(:lenpre)//'.CSA.bank'
1579 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1580 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1581 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1582 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1583 csa_int=prefix(:lenpre)//'.int'
1584 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1585 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1586 csa_in=prefix(:lenpre)//'.CSA.in'
1587 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1590 write (iout,'(80(1h-))')
1591 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1592 write (iout,'(80(1h-))')
1593 write (iout,*) "Input file : ",
1594 & pref_orig(:ilen(pref_orig))//'.inp'
1595 write (iout,*) "Output file : ",
1596 & outname(:ilen(outname))
1598 write (iout,*) "Sidechain potential file : ",
1599 & sidename(:ilen(sidename))
1601 write (iout,*) "SCp potential file : ",
1602 & scpname(:ilen(scpname))
1604 write (iout,*) "Electrostatic potential file : ",
1605 & elename(:ilen(elename))
1606 write (iout,*) "Cumulant coefficient file : ",
1607 & fouriername(:ilen(fouriername))
1608 write (iout,*) "Torsional parameter file : ",
1609 & torname(:ilen(torname))
1610 write (iout,*) "Double torsional parameter file : ",
1611 & tordname(:ilen(tordname))
1612 write (iout,*) "SCCOR parameter file : ",
1613 & sccorname(:ilen(sccorname))
1614 write (iout,*) "Bond & inertia constant file : ",
1615 & bondname(:ilen(bondname))
1616 write (iout,*) "Bending parameter file : ",
1617 & thetname(:ilen(thetname))
1618 write (iout,*) "Rotamer parameter file : ",
1619 & rotname(:ilen(rotname))
1620 write (iout,*) "Threading database : ",
1621 & patname(:ilen(patname))
1623 &write (iout,*)" DIRTMP : ",
1625 write (iout,'(80(1h-))')
1629 c----------------------------------------------------------------------------
1630 subroutine card_concat(card)
1631 implicit real*8 (a-h,o-z)
1632 include 'DIMENSIONS'
1633 include 'COMMON.IOUNITS'
1635 character*80 karta,ucase
1637 read (inp,'(a)') karta
1640 do while (karta(80:80).eq.'&')
1641 card=card(:ilen(card)+1)//karta(:79)
1642 read (inp,'(a)') karta
1645 card=card(:ilen(card)+1)//karta
1649 subroutine read_dist_constr
1650 implicit real*8 (a-h,o-z)
1651 include 'DIMENSIONS'
1655 include 'COMMON.SETUP'
1656 include 'COMMON.CONTROL'
1657 include 'COMMON.CHAIN'
1658 include 'COMMON.IOUNITS'
1659 include 'COMMON.SBRIDGE'
1660 integer ifrag_(2,100),ipair_(2,100)
1661 double precision wfrag_(100),wpair_(100)
1662 character*500 controlcard
1663 write (iout,*) "Calling read_dist_constr"
1664 write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1666 call card_concat(controlcard)
1667 call readi(controlcard,"NFRAG",nfrag_,0)
1668 call readi(controlcard,"NPAIR",npair_,0)
1669 call readi(controlcard,"NDIST",ndist_,0)
1670 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1671 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1672 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1673 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1674 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1675 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1676 c write (iout,*) "IFRAG"
1678 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1680 c write (iout,*) "IPAIR"
1682 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1686 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1687 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1688 & ifrag_(2,i)=nstart_sup+nsup-1
1689 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1691 if (wfrag_(i).gt.0.0d0) then
1692 do j=ifrag_(1,i),ifrag_(2,i)-1
1693 do k=j+1,ifrag_(2,i)
1694 write (iout,*) "j",j," k",k
1696 if (constr_dist.eq.1) then
1701 forcon(nhpb)=wfrag_(i)
1702 else if (constr_dist.eq.2) then
1703 if (ddjk.le.dist_cut) then
1708 forcon(nhpb)=wfrag_(i)
1715 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1718 if (.not.out1file .or. me.eq.king)
1719 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1720 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1722 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1723 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1730 if (wpair_(i).gt.0.0d0) then
1738 do j=ifrag_(1,ii),ifrag_(2,ii)
1739 do k=ifrag_(1,jj),ifrag_(2,jj)
1743 forcon(nhpb)=wpair_(i)
1744 dhpb(nhpb)=dist(j,k)
1746 if (.not.out1file .or. me.eq.king)
1747 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1748 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1750 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1751 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1758 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1759 if (forcon(nhpb+1).gt.0.0d0) then
1761 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1763 if (.not.out1file .or. me.eq.king)
1764 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1765 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1767 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1768 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1775 c-------------------------------------------------------------------------------
1777 subroutine flush(iu)
1782 subroutine flush(iu)
1787 c------------------------------------------------------------------------------
1788 subroutine copy_to_tmp(source)
1789 include "DIMENSIONS"
1790 include "COMMON.IOUNITS"
1791 character*(*) source
1792 character* 256 tmpfile
1796 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1797 inquire(file=tmpfile,exist=ex)
1799 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1800 & " to temporary directory..."
1801 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1802 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1806 c------------------------------------------------------------------------------
1807 subroutine move_from_tmp(source)
1808 include "DIMENSIONS"
1809 include "COMMON.IOUNITS"
1810 character*(*) source
1813 write (*,*) "Moving ",source(:ilen(source)),
1814 & " from temporary directory to working directory"
1815 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1816 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1819 c------------------------------------------------------------------------------
1820 subroutine random_init(seed)
1822 C Initialize random number generator
1824 implicit real*8 (a-h,o-z)
1825 include 'DIMENSIONS'
1831 logical OKRandom, prng_restart
1833 integer iseed_array(4)
1835 include 'COMMON.IOUNITS'
1836 include 'COMMON.TIME1'
1837 c include 'COMMON.THREAD'
1838 include 'COMMON.SBRIDGE'
1839 include 'COMMON.CONTROL'
1840 include 'COMMON.MCM'
1841 c include 'COMMON.MAP'
1842 include 'COMMON.HEADER'
1843 c include 'COMMON.CSA'
1844 include 'COMMON.CHAIN'
1845 c include 'COMMON.MUCA'
1846 c include 'COMMON.MD'
1847 include 'COMMON.FFIELD'
1848 include 'COMMON.SETUP'
1849 iseed=-dint(dabs(seed))
1850 if (iseed.eq.0) then
1851 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1852 & 'Random seed undefined. The program will stop.'
1853 write (*,'(/80(1h*)/20x,a/80(1h*))')
1854 & 'Random seed undefined. The program will stop.'
1856 call mpi_finalize(mpi_comm_world,ierr)
1858 stop 'Bad random seed.'
1861 if (fg_rank.eq.0) then
1865 if(me.eq.king .or. .not. out1file)
1866 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1867 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1868 OKRandom = prng_restart(me,iseedi8)
1871 tmp=65536.0d0**(4-i)
1872 iseed_array(i) = dint(seed/tmp)
1873 seed=seed-iseed_array(i)*tmp
1875 if(me.eq.king .or. .not. out1file)
1876 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1877 & (iseed_array(i),i=1,4)
1878 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1879 & (iseed_array(i),i=1,4)
1880 OKRandom = prng_restart(me,iseed_array)
1883 r1=ran_number(0.0D0,1.0D0)
1884 if(me.eq.king .or. .not. out1file)
1885 & write (iout,*) 'ran_num',r1
1886 if (r1.lt.0.0d0) OKRandom=.false.
1888 if (.not.OKRandom) then
1889 write (iout,*) 'PRNG IS NOT WORKING!!!'
1890 print *,'PRNG IS NOT WORKING!!!'
1893 call mpi_abort(mpi_comm_world,error_msg,ierr)
1896 write (iout,*) 'too many processors for parallel prng'
1897 write (*,*) 'too many processors for parallel prng'
1905 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)