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,0)
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 c 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 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1158 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1159 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1160 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1161 & 'Options in energy minimization:'
1162 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1163 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1164 & 'MinMin:',MinMin,' MinFun:',MinFun,
1165 & ' TolF:',TolF,' RTolF:',RTolF
1168 c----------------------------------------------------------------------------
1169 subroutine read_angles(kanal,*)
1170 implicit real*8 (a-h,o-z)
1171 include 'DIMENSIONS'
1172 include 'COMMON.GEO'
1173 include 'COMMON.VAR'
1174 include 'COMMON.CHAIN'
1175 include 'COMMON.IOUNITS'
1176 include 'COMMON.CONTROL'
1177 c Read angles from input
1179 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1180 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1181 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1182 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1185 c 9/7/01 avoid 180 deg valence angle
1186 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1188 theta(i)=deg2rad*theta(i)
1189 phi(i)=deg2rad*phi(i)
1190 alph(i)=deg2rad*alph(i)
1191 omeg(i)=deg2rad*omeg(i)
1196 c----------------------------------------------------------------------------
1197 subroutine reada(rekord,lancuch,wartosc,default)
1199 character*(*) rekord,lancuch
1200 double precision wartosc,default
1203 iread=index(rekord,lancuch)
1204 if (iread.eq.0) then
1208 iread=iread+ilen(lancuch)+1
1209 read (rekord(iread:),*,err=10,end=10) wartosc
1214 c----------------------------------------------------------------------------
1215 subroutine readi(rekord,lancuch,wartosc,default)
1217 character*(*) rekord,lancuch
1218 integer wartosc,default
1221 iread=index(rekord,lancuch)
1222 if (iread.eq.0) then
1226 iread=iread+ilen(lancuch)+1
1227 read (rekord(iread:),*,err=10,end=10) wartosc
1232 c----------------------------------------------------------------------------
1233 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1236 integer tablica(dim),default
1237 character*(*) rekord,lancuch
1244 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1245 if (iread.eq.0) return
1246 iread=iread+ilen(lancuch)+1
1247 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1250 c----------------------------------------------------------------------------
1251 subroutine multreada(rekord,lancuch,tablica,dim,default)
1254 double precision tablica(dim),default
1255 character*(*) rekord,lancuch
1262 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1263 if (iread.eq.0) return
1264 iread=iread+ilen(lancuch)+1
1265 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1268 c----------------------------------------------------------------------------
1269 subroutine openunits
1270 implicit real*8 (a-h,o-z)
1271 include 'DIMENSIONS'
1274 character*16 form,nodename
1277 include 'COMMON.SETUP'
1278 include 'COMMON.IOUNITS'
1279 c include 'COMMON.MD'
1280 include 'COMMON.CONTROL'
1281 integer lenpre,lenpot,ilen,lentmp
1283 character*3 out1file_text,ucase
1286 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1287 call getenv_loc("PREFIX",prefix)
1289 call getenv_loc("POT",pot)
1290 call getenv_loc("DIRTMP",tmpdir)
1291 call getenv_loc("CURDIR",curdir)
1292 call getenv_loc("OUT1FILE",out1file_text)
1293 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1294 out1file_text=ucase(out1file_text)
1295 if (out1file_text(1:1).eq."Y") then
1298 out1file=fg_rank.gt.0
1303 if (lentmp.gt.0) then
1304 write (*,'(80(1h!))')
1305 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1306 write (*,'(80(1h!))')
1307 write (*,*)"All output files will be on node /tmp directory."
1309 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1310 if (me.eq.king) then
1311 write (*,*) "The master node is ",nodename
1312 else if (fg_rank.eq.0) then
1313 write (*,*) "I am the CG slave node ",nodename
1315 write (*,*) "I am the FG slave node ",nodename
1318 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1319 lenpre = lentmp+lenpre+1
1321 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1322 C Get the names and open the input files
1323 #if defined(WINIFL) || defined(WINPGI)
1324 open(1,file=pref_orig(:ilen(pref_orig))//
1325 & '.inp',status='old',readonly,shared)
1326 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1327 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1328 C Get parameter filenames and open the parameter files.
1329 call getenv_loc('BONDPAR',bondname)
1330 open (ibond,file=bondname,status='old',readonly,shared)
1331 call getenv_loc('THETPAR',thetname)
1332 open (ithep,file=thetname,status='old',readonly,shared)
1334 call getenv_loc('THETPARPDB',thetname_pdb)
1335 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1337 call getenv_loc('ROTPAR',rotname)
1338 open (irotam,file=rotname,status='old',readonly,shared)
1340 call getenv_loc('ROTPARPDB',rotname_pdb)
1341 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1343 call getenv_loc('TORPAR',torname)
1344 open (itorp,file=torname,status='old',readonly,shared)
1345 call getenv_loc('TORDPAR',tordname)
1346 open (itordp,file=tordname,status='old',readonly,shared)
1347 call getenv_loc('FOURIER',fouriername)
1348 open (ifourier,file=fouriername,status='old',readonly,shared)
1349 call getenv_loc('ELEPAR',elename)
1350 open (ielep,file=elename,status='old',readonly,shared)
1351 call getenv_loc('SIDEPAR',sidename)
1352 open (isidep,file=sidename,status='old',readonly,shared)
1353 #elif (defined CRAY) || (defined AIX)
1354 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1356 c print *,"Processor",myrank," opened file 1"
1357 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1358 c print *,"Processor",myrank," opened file 9"
1359 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1360 C Get parameter filenames and open the parameter files.
1361 call getenv_loc('BONDPAR',bondname)
1362 open (ibond,file=bondname,status='old',action='read')
1363 c print *,"Processor",myrank," opened file IBOND"
1364 call getenv_loc('THETPAR',thetname)
1365 open (ithep,file=thetname,status='old',action='read')
1366 c print *,"Processor",myrank," opened file ITHEP"
1368 call getenv_loc('THETPARPDB',thetname_pdb)
1369 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1371 call getenv_loc('ROTPAR',rotname)
1372 open (irotam,file=rotname,status='old',action='read')
1373 c print *,"Processor",myrank," opened file IROTAM"
1375 call getenv_loc('ROTPARPDB',rotname_pdb)
1376 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1378 call getenv_loc('TORPAR',torname)
1379 open (itorp,file=torname,status='old',action='read')
1380 c print *,"Processor",myrank," opened file ITORP"
1381 call getenv_loc('TORDPAR',tordname)
1382 open (itordp,file=tordname,status='old',action='read')
1383 c print *,"Processor",myrank," opened file ITORDP"
1384 call getenv_loc('SCCORPAR',sccorname)
1385 open (isccor,file=sccorname,status='old',action='read')
1386 c print *,"Processor",myrank," opened file ISCCOR"
1387 call getenv_loc('FOURIER',fouriername)
1388 open (ifourier,file=fouriername,status='old',action='read')
1389 c print *,"Processor",myrank," opened file IFOURIER"
1390 call getenv_loc('ELEPAR',elename)
1391 open (ielep,file=elename,status='old',action='read')
1392 c print *,"Processor",myrank," opened file IELEP"
1393 call getenv_loc('SIDEPAR',sidename)
1394 open (isidep,file=sidename,status='old',action='read')
1395 c print *,"Processor",myrank," opened file ISIDEP"
1396 c print *,"Processor",myrank," opened parameter files"
1398 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1399 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1400 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1401 C Get parameter filenames and open the parameter files.
1402 call getenv_loc('BONDPAR',bondname)
1403 open (ibond,file=bondname,status='old')
1404 call getenv_loc('THETPAR',thetname)
1405 open (ithep,file=thetname,status='old')
1407 call getenv_loc('THETPARPDB',thetname_pdb)
1408 open (ithep_pdb,file=thetname_pdb,status='old')
1410 call getenv_loc('ROTPAR',rotname)
1411 open (irotam,file=rotname,status='old')
1413 call getenv_loc('ROTPARPDB',rotname_pdb)
1414 open (irotam_pdb,file=rotname_pdb,status='old')
1416 call getenv_loc('TORPAR',torname)
1417 open (itorp,file=torname,status='old')
1418 call getenv_loc('TORDPAR',tordname)
1419 open (itordp,file=tordname,status='old')
1420 call getenv_loc('SCCORPAR',sccorname)
1421 open (isccor,file=sccorname,status='old')
1422 call getenv_loc('FOURIER',fouriername)
1423 open (ifourier,file=fouriername,status='old')
1424 call getenv_loc('ELEPAR',elename)
1425 open (ielep,file=elename,status='old')
1426 call getenv_loc('SIDEPAR',sidename)
1427 open (isidep,file=sidename,status='old')
1429 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1431 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1432 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1433 C Get parameter filenames and open the parameter files.
1434 call getenv_loc('BONDPAR',bondname)
1435 open (ibond,file=bondname,status='old',readonly)
1436 call getenv_loc('THETPAR',thetname)
1437 open (ithep,file=thetname,status='old',readonly)
1439 call getenv_loc('THETPARPDB',thetname_pdb)
1440 print *,"thetname_pdb ",thetname_pdb
1441 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1442 print *,ithep_pdb," opened"
1444 call getenv_loc('ROTPAR',rotname)
1445 open (irotam,file=rotname,status='old',readonly)
1447 call getenv_loc('ROTPARPDB',rotname_pdb)
1448 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1450 call getenv_loc('TORPAR',torname)
1451 open (itorp,file=torname,status='old',readonly)
1452 call getenv_loc('TORDPAR',tordname)
1453 open (itordp,file=tordname,status='old',readonly)
1454 call getenv_loc('SCCORPAR',sccorname)
1455 open (isccor,file=sccorname,status='old',readonly)
1456 call getenv_loc('FOURIER',fouriername)
1457 open (ifourier,file=fouriername,status='old',readonly)
1458 call getenv_loc('ELEPAR',elename)
1459 open (ielep,file=elename,status='old',readonly)
1460 call getenv_loc('SIDEPAR',sidename)
1461 open (isidep,file=sidename,status='old',readonly)
1465 C 8/9/01 In the newest version SCp interaction constants are read from a file
1466 C Use -DOLDSCP to use hard-coded constants instead.
1468 call getenv_loc('SCPPAR',scpname)
1469 #if defined(WINIFL) || defined(WINPGI)
1470 open (iscpp,file=scpname,status='old',readonly,shared)
1471 #elif (defined CRAY) || (defined AIX)
1472 open (iscpp,file=scpname,status='old',action='read')
1474 open (iscpp,file=scpname,status='old')
1476 open (iscpp,file=scpname,status='old',readonly)
1479 call getenv_loc('PATTERN',patname)
1480 #if defined(WINIFL) || defined(WINPGI)
1481 open (icbase,file=patname,status='old',readonly,shared)
1482 #elif (defined CRAY) || (defined AIX)
1483 open (icbase,file=patname,status='old',action='read')
1485 open (icbase,file=patname,status='old')
1487 open (icbase,file=patname,status='old',readonly)
1490 C Open output file only for CG processes
1491 c print *,"Processor",myrank," fg_rank",fg_rank
1492 if (fg_rank.eq.0) then
1494 if (nodes.eq.1) then
1497 npos = dlog10(dfloat(nodes-1))+1
1499 if (npos.lt.3) npos=3
1500 write (liczba,'(i1)') npos
1501 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1503 write (liczba,form) me
1504 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1505 & liczba(:ilen(liczba))
1506 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1508 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1510 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1511 & liczba(:ilen(liczba))//'.mol2'
1512 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1513 & liczba(:ilen(liczba))//'.stat'
1515 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1516 & //liczba(:ilen(liczba))//'.stat')
1517 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1520 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1521 c & liczba(:ilen(liczba))//'.const'
1526 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1527 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1528 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1529 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1530 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1532 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1534 rest2name=prefix(:ilen(prefix))//'.rst'
1536 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1539 #if defined(AIX) || defined(PGI)
1540 if (me.eq.king .or. .not. out1file)
1541 & open(iout,file=outname,status='unknown')
1543 if (fg_rank.gt.0) then
1544 write (liczba,'(i3.3)') myrank/nfgtasks
1545 write (ll,'(bz,i3.3)') fg_rank
1546 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1551 open(igeom,file=intname,status='unknown',position='append')
1552 open(ipdb,file=pdbname,status='unknown')
1553 open(imol2,file=mol2name,status='unknown')
1554 open(istat,file=statname,status='unknown',position='append')
1556 c1out open(iout,file=outname,status='unknown')
1559 if (me.eq.king .or. .not.out1file)
1560 & open(iout,file=outname,status='unknown')
1562 if (fg_rank.gt.0) then
1563 write (liczba,'(i3.3)') myrank/nfgtasks
1564 write (ll,'(bz,i3.3)') fg_rank
1565 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1570 open(igeom,file=intname,status='unknown',access='append')
1571 open(ipdb,file=pdbname,status='unknown')
1572 open(imol2,file=mol2name,status='unknown')
1573 open(istat,file=statname,status='unknown',access='append')
1575 c1out open(iout,file=outname,status='unknown')
1578 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1579 csa_seed=prefix(:lenpre)//'.CSA.seed'
1580 csa_history=prefix(:lenpre)//'.CSA.history'
1581 csa_bank=prefix(:lenpre)//'.CSA.bank'
1582 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1583 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1584 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1585 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1586 csa_int=prefix(:lenpre)//'.int'
1587 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1588 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1589 csa_in=prefix(:lenpre)//'.CSA.in'
1590 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1593 write (iout,'(80(1h-))')
1594 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1595 write (iout,'(80(1h-))')
1596 write (iout,*) "Input file : ",
1597 & pref_orig(:ilen(pref_orig))//'.inp'
1598 write (iout,*) "Output file : ",
1599 & outname(:ilen(outname))
1601 write (iout,*) "Sidechain potential file : ",
1602 & sidename(:ilen(sidename))
1604 write (iout,*) "SCp potential file : ",
1605 & scpname(:ilen(scpname))
1607 write (iout,*) "Electrostatic potential file : ",
1608 & elename(:ilen(elename))
1609 write (iout,*) "Cumulant coefficient file : ",
1610 & fouriername(:ilen(fouriername))
1611 write (iout,*) "Torsional parameter file : ",
1612 & torname(:ilen(torname))
1613 write (iout,*) "Double torsional parameter file : ",
1614 & tordname(:ilen(tordname))
1615 write (iout,*) "SCCOR parameter file : ",
1616 & sccorname(:ilen(sccorname))
1617 write (iout,*) "Bond & inertia constant file : ",
1618 & bondname(:ilen(bondname))
1619 write (iout,*) "Bending parameter file : ",
1620 & thetname(:ilen(thetname))
1621 write (iout,*) "Rotamer parameter file : ",
1622 & rotname(:ilen(rotname))
1623 write (iout,*) "Threading database : ",
1624 & patname(:ilen(patname))
1626 &write (iout,*)" DIRTMP : ",
1628 write (iout,'(80(1h-))')
1632 c----------------------------------------------------------------------------
1633 subroutine card_concat(card)
1634 implicit real*8 (a-h,o-z)
1635 include 'DIMENSIONS'
1636 include 'COMMON.IOUNITS'
1638 character*80 karta,ucase
1640 read (inp,'(a)') karta
1643 do while (karta(80:80).eq.'&')
1644 card=card(:ilen(card)+1)//karta(:79)
1645 read (inp,'(a)') karta
1648 card=card(:ilen(card)+1)//karta
1652 subroutine read_dist_constr
1653 implicit real*8 (a-h,o-z)
1654 include 'DIMENSIONS'
1658 include 'COMMON.SETUP'
1659 include 'COMMON.CONTROL'
1660 include 'COMMON.CHAIN'
1661 include 'COMMON.IOUNITS'
1662 include 'COMMON.SBRIDGE'
1663 integer ifrag_(2,100),ipair_(2,100)
1664 double precision wfrag_(100),wpair_(100)
1665 character*500 controlcard
1666 write (iout,*) "Calling read_dist_constr"
1667 write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1669 call card_concat(controlcard)
1670 call readi(controlcard,"NFRAG",nfrag_,0)
1671 call readi(controlcard,"NPAIR",npair_,0)
1672 call readi(controlcard,"NDIST",ndist_,0)
1673 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1674 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1675 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1676 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1677 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1678 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1679 c write (iout,*) "IFRAG"
1681 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1683 c write (iout,*) "IPAIR"
1685 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1689 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1690 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1691 & ifrag_(2,i)=nstart_sup+nsup-1
1692 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1694 if (wfrag_(i).gt.0.0d0) then
1695 do j=ifrag_(1,i),ifrag_(2,i)-1
1696 do k=j+1,ifrag_(2,i)
1697 write (iout,*) "j",j," k",k
1699 if (constr_dist.eq.1) then
1704 forcon(nhpb)=wfrag_(i)
1705 else if (constr_dist.eq.2) then
1706 if (ddjk.le.dist_cut) then
1711 forcon(nhpb)=wfrag_(i)
1718 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1721 if (.not.out1file .or. me.eq.king)
1722 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1723 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1725 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1726 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1733 if (wpair_(i).gt.0.0d0) then
1741 do j=ifrag_(1,ii),ifrag_(2,ii)
1742 do k=ifrag_(1,jj),ifrag_(2,jj)
1746 forcon(nhpb)=wpair_(i)
1747 dhpb(nhpb)=dist(j,k)
1749 if (.not.out1file .or. me.eq.king)
1750 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1751 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1753 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1754 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1761 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1762 if (forcon(nhpb+1).gt.0.0d0) then
1764 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1766 if (.not.out1file .or. me.eq.king)
1767 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1768 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1770 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1771 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1778 c-------------------------------------------------------------------------------
1780 subroutine flush(iu)
1785 subroutine flush(iu)
1790 c------------------------------------------------------------------------------
1791 subroutine copy_to_tmp(source)
1792 include "DIMENSIONS"
1793 include "COMMON.IOUNITS"
1794 character*(*) source
1795 character* 256 tmpfile
1799 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1800 inquire(file=tmpfile,exist=ex)
1802 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1803 & " to temporary directory..."
1804 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1805 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1809 c------------------------------------------------------------------------------
1810 subroutine move_from_tmp(source)
1811 include "DIMENSIONS"
1812 include "COMMON.IOUNITS"
1813 character*(*) source
1816 write (*,*) "Moving ",source(:ilen(source)),
1817 & " from temporary directory to working directory"
1818 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1819 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1822 c------------------------------------------------------------------------------
1823 subroutine random_init(seed)
1825 C Initialize random number generator
1827 implicit real*8 (a-h,o-z)
1828 include 'DIMENSIONS'
1834 logical OKRandom, prng_restart
1836 integer iseed_array(4)
1838 include 'COMMON.IOUNITS'
1839 include 'COMMON.TIME1'
1840 c include 'COMMON.THREAD'
1841 include 'COMMON.SBRIDGE'
1842 include 'COMMON.CONTROL'
1843 include 'COMMON.MCM'
1844 c include 'COMMON.MAP'
1845 include 'COMMON.HEADER'
1846 c include 'COMMON.CSA'
1847 include 'COMMON.CHAIN'
1848 c include 'COMMON.MUCA'
1849 c include 'COMMON.MD'
1850 include 'COMMON.FFIELD'
1851 include 'COMMON.SETUP'
1852 iseed=-dint(dabs(seed))
1853 if (iseed.eq.0) then
1854 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1855 & 'Random seed undefined. The program will stop.'
1856 write (*,'(/80(1h*)/20x,a/80(1h*))')
1857 & 'Random seed undefined. The program will stop.'
1859 call mpi_finalize(mpi_comm_world,ierr)
1861 stop 'Bad random seed.'
1864 if (fg_rank.eq.0) then
1868 if(me.eq.king .or. .not. out1file)
1869 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1870 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1871 OKRandom = prng_restart(me,iseedi8)
1874 tmp=65536.0d0**(4-i)
1875 iseed_array(i) = dint(seed/tmp)
1876 seed=seed-iseed_array(i)*tmp
1878 if(me.eq.king .or. .not. out1file)
1879 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1880 & (iseed_array(i),i=1,4)
1881 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1882 & (iseed_array(i),i=1,4)
1883 OKRandom = prng_restart(me,iseed_array)
1886 r1=ran_number(0.0D0,1.0D0)
1887 if(me.eq.king .or. .not. out1file)
1888 & write (iout,*) 'ran_num',r1
1889 if (r1.lt.0.0d0) OKRandom=.false.
1891 if (.not.OKRandom) then
1892 write (iout,*) 'PRNG IS NOT WORKING!!!'
1893 print *,'PRNG IS NOT WORKING!!!'
1896 call mpi_abort(mpi_comm_world,error_msg,ierr)
1899 write (iout,*) 'too many processors for parallel prng'
1900 write (*,*) 'too many processors for parallel prng'
1908 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)