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'
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
91 read (INP,'(a)') titel
92 call card_concat(controlcard)
93 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
94 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
95 call reada(controlcard,'SEED',seed,0.0D0)
96 call random_init(seed)
97 C Set up the time limit (caution! The time must be input in minutes!)
98 read_cart=index(controlcard,'READ_CART').gt.0
99 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
200 lmuca=index(controlcard,'MUCA').gt.0
201 call readi(controlcard,'MUCADYN',mucadyn,0)
202 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203 if (lmuca .and. (me.eq.king .or. .not.out1file ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
209 iscode=index(controlcard,'ONE_LETTER')
210 indphi=index(controlcard,'PHI')
211 indback=index(controlcard,'BACK')
212 iranconf=index(controlcard,'RAND_CONF')
213 i2ndstr=index(controlcard,'USE_SEC_PRED')
214 gradout=index(controlcard,'GRADOUT').gt.0
215 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
217 if(me.eq.king.or..not.out1file)
218 & write (iout,'(2a)') diagmeth(kdiag),
219 & ' routine used to diagonalize matrices.'
222 c------------------------------------------------------------------------------
225 C Read molecular data.
227 implicit real*8 (a-h,o-z)
233 include 'COMMON.IOUNITS'
236 include 'COMMON.INTERACT'
237 include 'COMMON.LOCAL'
238 include 'COMMON.NAMES'
239 include 'COMMON.CHAIN'
240 include 'COMMON.FFIELD'
241 include 'COMMON.SBRIDGE'
242 include 'COMMON.HEADER'
243 include 'COMMON.CONTROL'
244 c include 'COMMON.DBASE'
245 c include 'COMMON.THREAD'
246 include 'COMMON.CONTACTS'
247 include 'COMMON.TORCNSTR'
248 include 'COMMON.TIME1'
249 include 'COMMON.BOUNDS'
250 c include 'COMMON.MD'
251 c include 'COMMON.REMD'
252 include 'COMMON.SETUP'
253 character*4 sequence(maxres)
255 double precision x(maxvar)
256 character*256 pdbfile
257 character*320 weightcard
258 character*80 weightcard_t,ucase
259 dimension itype_pdb(maxres)
260 common /pizda/ itype_pdb
261 logical seq_comp,fail
262 double precision energia(0:n_ene)
268 C Read weights of the subsequent energy terms.
270 call card_concat(weightcard)
271 call reada(weightcard,'WLONG',wlong,1.0D0)
272 call reada(weightcard,'WSC',wsc,wlong)
273 call reada(weightcard,'WSCP',wscp,wlong)
274 call reada(weightcard,'WELEC',welec,1.0D0)
275 call reada(weightcard,'WVDWPP',wvdwpp,welec)
276 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
277 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
278 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
279 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
280 call reada(weightcard,'WTURN3',wturn3,1.0D0)
281 call reada(weightcard,'WTURN4',wturn4,1.0D0)
282 call reada(weightcard,'WTURN6',wturn6,1.0D0)
283 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
284 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
285 call reada(weightcard,'WBOND',wbond,1.0D0)
286 call reada(weightcard,'WTOR',wtor,1.0D0)
287 call reada(weightcard,'WTORD',wtor_d,1.0D0)
288 call reada(weightcard,'WANG',wang,1.0D0)
289 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
291 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
292 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
293 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
294 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
296 call reada(weightcard,'SCAL14',scal14,0.4D0)
297 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
298 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
299 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
300 call reada(weightcard,'TEMP0',temp0,300.0d0)
301 if (index(weightcard,'SOFT').gt.0) ipot=6
302 C 12/1/95 Added weight for the multi-body term WCORR
303 call reada(weightcard,'WCORRH',wcorr,1.0D0)
304 if (wcorr4.gt.0.0d0) wcorr=wcorr4
325 weights(24)=wdfa_dist
328 weights(27)=wdfa_beta
331 if(me.eq.king.or..not.out1file)
332 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
333 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
335 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
337 10 format (/'Energy-term weights (unscaled):'//
338 & 'WSCC= ',f10.6,' (SC-SC)'/
339 & 'WSCP= ',f10.6,' (SC-p)'/
340 & 'WELEC= ',f10.6,' (p-p electr)'/
341 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
342 & 'WBOND= ',f10.6,' (stretching)'/
343 & 'WANG= ',f10.6,' (bending)'/
344 & 'WSCLOC= ',f10.6,' (SC local)'/
345 & 'WTOR= ',f10.6,' (torsional)'/
346 & 'WTORD= ',f10.6,' (double torsional)'/
347 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
348 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
349 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
350 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
351 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
352 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
353 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
354 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
355 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
356 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
357 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
358 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
359 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
361 if(me.eq.king.or..not.out1file)then
362 if (wcorr4.gt.0.0d0) then
363 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
364 & 'between contact pairs of peptide groups'
365 write (iout,'(2(a,f5.3/))')
366 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
367 & 'Range of quenching the correlation terms:',2*delt_corr
368 else if (wcorr.gt.0.0d0) then
369 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
370 & 'between contact pairs of peptide groups'
372 write (iout,'(a,f8.3)')
373 & 'Scaling factor of 1,4 SC-p interactions:',scal14
374 write (iout,'(a,f8.3)')
375 & 'General scaling factor of SC-p interactions:',scalscp
377 r0_corr=cutoff_corr-delt_corr
379 aad(i,1)=scalscp*aad(i,1)
380 aad(i,2)=scalscp*aad(i,2)
381 bad(i,1)=scalscp*bad(i,1)
382 bad(i,2)=scalscp*bad(i,2)
384 c call rescale_weights(t_bath)
385 if(me.eq.king.or..not.out1file)
386 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
387 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
389 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
391 22 format (/'Energy-term weights (scaled):'//
392 & 'WSCC= ',f10.6,' (SC-SC)'/
393 & 'WSCP= ',f10.6,' (SC-p)'/
394 & 'WELEC= ',f10.6,' (p-p electr)'/
395 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
396 & 'WBOND= ',f10.6,' (stretching)'/
397 & 'WANG= ',f10.6,' (bending)'/
398 & 'WSCLOC= ',f10.6,' (SC local)'/
399 & 'WTOR= ',f10.6,' (torsional)'/
400 & 'WTORD= ',f10.6,' (double torsional)'/
401 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
402 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
403 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
404 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
405 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
406 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
407 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
408 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
409 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
410 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
411 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
412 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
413 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
415 if(me.eq.king.or..not.out1file)
416 & write (iout,*) "Reference temperature for weights calculation:",
418 call reada(weightcard,"D0CM",d0cm,3.78d0)
419 call reada(weightcard,"AKCM",akcm,15.1d0)
420 call reada(weightcard,"AKTH",akth,11.0d0)
421 call reada(weightcard,"AKCT",akct,12.0d0)
422 call reada(weightcard,"V1SS",v1ss,-1.08d0)
423 call reada(weightcard,"V2SS",v2ss,7.61d0)
424 call reada(weightcard,"V3SS",v3ss,13.7d0)
425 call reada(weightcard,"EBR",ebr,-5.50D0)
426 if(me.eq.king.or..not.out1file) then
427 write (iout,*) "Parameters of the SS-bond potential:"
428 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
430 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
431 write (iout,*) "EBR",ebr
432 print *,'indpdb=',indpdb,' pdbref=',pdbref
434 if (indpdb.gt.0 .or. pdbref) then
435 read(inp,'(a)') pdbfile
436 if(me.eq.king.or..not.out1file)
437 & write (iout,'(2a)') 'PDB data will be read from file ',
438 & pdbfile(:ilen(pdbfile))
439 open(ipdbin,file=pdbfile,status='old',err=33)
441 33 write (iout,'(a)') 'Error opening PDB file.'
444 c print *,'Begin reading pdb data'
446 c print *,'Finished reading pdb data'
447 if(me.eq.king.or..not.out1file)
448 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
449 & ' nstart_sup=',nstart_sup
451 itype_pdb(i)=itype(i)
455 nct=nstart_sup+nsup-1
456 call contact(.false.,ncont_ref,icont_ref,co)
459 C Following 2 lines for diagnostics; comment out if not needed
460 write (iout,*) "Before sideadd"
462 if(me.eq.king.or..not.out1file)
463 & write(iout,*)'Adding sidechains'
470 do while (fail.and.nsi.le.maxsi)
471 c call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
474 if(fail) write(iout,*)'Adding sidechain failed for res ',
475 & i,' after ',nsi,' trials'
479 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
481 C Following 2 lines for diagnostics; comment out if not needed
482 write (iout,*) "After sideadd"
486 if (indpdb.eq.0) then
487 C Read sequence if not taken from the pdb file.
489 c print *,'nres=',nres
490 if (iscode.gt.0) then
491 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
493 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
495 C Convert sequence to numeric code
497 itype(i)=rescode(i,sequence(i),iscode)
499 C Assign initial virtual bond lengths
505 vbld(i+nres)=dsc(itype(i))
506 vbld_inv(i+nres)=dsc_inv(itype(i))
507 c write (iout,*) "i",i," itype",itype(i),
508 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
512 c print '(20i4)',(itype(i),i=1,nres)
515 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
517 if (itype(i).eq.21) then
521 else if (itype(i+1).ne.20) then
523 else if (itype(i).ne.20) then
530 if(me.eq.king.or..not.out1file)then
531 write (iout,*) "ITEL"
533 write (iout,*) i,itype(i),itel(i)
535 print *,'Call Read_Bridge.'
538 C 8/13/98 Set limits to generating the dihedral angles
543 read (inp,*) ndih_constr
544 if (ndih_constr.gt.0) then
546 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
547 if(me.eq.king.or..not.out1file)then
549 & 'There are',ndih_constr,' constraints on phi angles.'
551 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
555 phi0(i)=deg2rad*phi0(i)
556 drange(i)=deg2rad*drange(i)
558 if(me.eq.king.or..not.out1file)
559 & write (iout,*) 'FTORS',ftors
562 phibound(1,ii) = phi0(i)-drange(i)
563 phibound(2,ii) = phi0(i)+drange(i)
570 write (iout,'(a)') 'Boundaries in phi angle sampling:'
572 write (iout,'(a3,i5,2f10.1)')
573 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
579 cd print *,'NNT=',NNT,' NCT=',NCT
580 if (itype(1).eq.21) nnt=2
581 if (itype(nres).eq.21) nct=nct-1
583 C Juyong:READ init_vars
584 C Initialize variables!
585 C Juyong:READ read_info
586 C READ fragment information!!
587 C both routines should be in dfa.F file!!
589 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
590 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
592 print*, 'init_dfa_vars finished!'
594 print*, 'read_dfa_info finished!'
601 if(me.eq.king.or..not.out1file)
602 & write (iout,'(a,i3)') 'nsup=',nsup
604 if (nsup.le.(nct-nnt+1)) then
605 do i=0,nct-nnt+1-nsup
606 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
612 & 'Error - sequences to be superposed do not match.'
615 do i=0,nsup-(nct-nnt+1)
616 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
618 nstart_sup=nstart_sup+i
624 & 'Error - sequences to be superposed do not match.'
627 if (nsup.eq.0) nsup=nct-nnt
628 if (nstart_sup.eq.0) nstart_sup=nnt
629 if (nstart_seq.eq.0) nstart_seq=nnt
630 if(me.eq.king.or..not.out1file)
631 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
632 & ' nstart_seq=',nstart_seq
634 c--- Zscore rms -------
635 if (nz_start.eq.0) nz_start=nnt
636 if (nz_end.eq.0 .and. nsup.gt.0) then
638 else if (nz_end.eq.0) then
641 if(me.eq.king.or..not.out1file)then
642 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
643 write (iout,*) 'IZ_SC=',iz_sc
645 c----------------------
648 if (.not.pdbref) then
649 call read_angles(inp,*38)
651 38 write (iout,'(a)') 'Error reading reference structure.'
653 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
654 stop 'Error reading reference structure'
658 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
667 call contact(.true.,ncont_ref,icont_ref,co)
669 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
671 if (constr_dist.gt.0) call read_dist_constr
672 c write (iout,*) "After read_dist_constr nhpb",nhpb
674 if(me.eq.king.or..not.out1file)
675 & write (iout,*) 'Contact order:',co
677 if(me.eq.king.or..not.out1file)
678 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
681 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
683 if(me.eq.king.or..not.out1file)
684 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
685 & icont_ref(1,i),' ',
686 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
690 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
691 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
692 & modecalc.ne.10) then
693 C If input structure hasn't been supplied from the PDB file read or generate
695 if (iranconf.eq.0 .and. .not. extconf) then
696 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
697 & write (iout,'(a)') 'Initial geometry will be read in.'
699 read(inp,'(8f10.5)',end=36,err=36)
700 & ((c(l,k),l=1,3),k=1,nres),
701 & ((c(l,k+nres),l=1,3),k=nnt,nct)
702 call int_from_cart1(.false.)
705 dc(j,i)=c(j,i+1)-c(j,i)
706 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
710 if (itype(i).ne.10) then
712 dc(j,i+nres)=c(j,i+nres)-c(j,i)
713 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
719 call read_angles(inp,*36)
722 36 write (iout,'(a)') 'Error reading angle file.'
724 call mpi_finalize( MPI_COMM_WORLD,IERR )
726 stop 'Error reading angle file.'
728 else if (extconf) then
729 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
730 & write (iout,'(a)') 'Extended chain initial geometry.'
732 theta(i)=90d0*deg2rad
738 alph(i)=110d0*deg2rad
741 omeg(i)=-120d0*deg2rad
744 if(me.eq.king.or..not.out1file)
745 & write (iout,'(a)') 'Random-generated initial geometry.'
749 if (me.eq.king .or. fg_rank.eq.0 .and. (
750 & modecalc.eq.12 .or. modecalc.eq.14) ) then
754 c call gen_rand_conf(itmp,*30)
756 30 write (iout,*) 'Failed to generate random conformation',
758 write (*,*) 'Processor:',me,
759 & ' Failed to generate random conformation',
769 write (iout,'(a,i3,a)') 'Processor:',me,
770 & ' error in generating random conformation.'
771 write (*,'(a,i3,a)') 'Processor:',me,
772 & ' error in generating random conformation.'
775 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
781 c call gen_rand_conf(itmp,*31)
783 31 write (iout,*) 'Failed to generate random conformation',
785 write (*,*) 'Failed to generate random conformation',
788 write (iout,'(a,i3,a)') 'Processor:',me,
789 & ' error in generating random conformation.'
790 write (*,'(a,i3,a)') 'Processor:',me,
791 & ' error in generating random conformation.'
796 elseif (modecalc.eq.4) then
797 read (inp,'(a)') intinname
798 open (intin,file=intinname,status='old',err=333)
799 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
800 & write (iout,'(a)') 'intinname',intinname
801 write (*,'(a)') 'Processor',myrank,' intinname',intinname
803 333 write (iout,'(2a)') 'Error opening angle file ',intinname
805 call MPI_Finalize(MPI_COMM_WORLD,IERR)
807 stop 'Error opening angle file.'
811 C Generate distance constraints, if the PDB structure is to be regularized.
813 if (me.eq.king .or. .not. out1file)
815 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
816 write (iout,'(/a,i3,a)')
817 & 'The chain contains',ns,' disulfide-bridging cysteines.'
818 write (iout,'(20i4)') (iss(i),i=1,ns)
819 write (iout,'(/a/)') 'Pre-formed links are:'
825 if (me.eq.king.or..not.out1file)
826 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
827 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
832 c if (i2ndstr.gt.0) call secstrp2dihc
833 c call geom_to_var(nvar,x)
834 c call etotal(energia(0))
835 c call enerprint(energia(0))
836 c call briefout(0,etot)
838 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
839 cd write (iout,'(a)') 'Variable list:'
840 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
842 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
843 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
844 & 'Processor',myrank,': end reading molecular data.'
848 c--------------------------------------------------------------------------
849 logical function seq_comp(itypea,itypeb,length)
851 integer length,itypea(length),itypeb(length)
854 if (itypea(i).ne.itypeb(i)) then
862 c-----------------------------------------------------------------------------
863 subroutine read_bridge
864 C Read information about disulfide bridges.
865 implicit real*8 (a-h,o-z)
870 include 'COMMON.IOUNITS'
873 include 'COMMON.INTERACT'
874 include 'COMMON.LOCAL'
875 include 'COMMON.NAMES'
876 include 'COMMON.CHAIN'
877 include 'COMMON.FFIELD'
878 include 'COMMON.SBRIDGE'
879 include 'COMMON.HEADER'
880 include 'COMMON.CONTROL'
881 c include 'COMMON.DBASE'
882 c include 'COMMON.THREAD'
883 include 'COMMON.TIME1'
884 include 'COMMON.SETUP'
885 C Read bridging residues.
886 read (inp,*) ns,(iss(i),i=1,ns)
888 if(me.eq.king.or..not.out1file)
889 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
890 C Check whether the specified bridging residues are cystines.
892 if (itype(iss(i)).ne.1) then
893 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
894 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
895 & ' can form a disulfide bridge?!!!'
896 write (*,'(2a,i3,a)')
897 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
898 & ' can form a disulfide bridge?!!!'
900 call MPI_Finalize(MPI_COMM_WORLD,ierror)
905 C Read preformed bridges.
907 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
908 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
911 C Check if the residues involved in bridges are in the specified list of
915 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
916 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
917 write (iout,'(a,i3,a)') 'Disulfide pair',i,
918 & ' contains residues present in other pairs.'
919 write (*,'(a,i3,a)') 'Disulfide pair',i,
920 & ' contains residues present in other pairs.'
922 call MPI_Finalize(MPI_COMM_WORLD,ierror)
928 if (ihpb(i).eq.iss(j)) goto 10
930 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
933 if (jhpb(i).eq.iss(j)) goto 20
935 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
948 c----------------------------------------------------------------------------
949 subroutine read_x(kanal,*)
950 implicit real*8 (a-h,o-z)
954 include 'COMMON.CHAIN'
955 include 'COMMON.IOUNITS'
956 include 'COMMON.CONTROL'
957 include 'COMMON.LOCAL'
958 include 'COMMON.INTERACT'
959 c Read coordinates from input
961 read(kanal,'(8f10.5)',end=10,err=10)
962 & ((c(l,k),l=1,3),k=1,nres),
963 & ((c(l,k+nres),l=1,3),k=nnt,nct)
966 c(j,2*nres)=c(j,nres)
968 call int_from_cart1(.false.)
971 dc(j,i)=c(j,i+1)-c(j,i)
972 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
976 if (itype(i).ne.10) then
978 dc(j,i+nres)=c(j,i+nres)-c(j,i)
979 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
987 c------------------------------------------------------------------------------
989 implicit real*8 (a-h,o-z)
991 include 'COMMON.IOUNITS'
994 include 'COMMON.INTERACT'
995 include 'COMMON.LOCAL'
996 include 'COMMON.NAMES'
997 include 'COMMON.CHAIN'
998 include 'COMMON.FFIELD'
999 include 'COMMON.SBRIDGE'
1000 include 'COMMON.HEADER'
1001 include 'COMMON.CONTROL'
1002 c include 'COMMON.DBASE'
1003 c include 'COMMON.THREAD'
1004 include 'COMMON.TIME1'
1005 C Set up variable list.
1011 if (itype(i).ne.10) then
1013 ialph(i,1)=nvar+nside
1017 if (indphi.gt.0) then
1019 else if (indback.gt.0) then
1024 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1027 c----------------------------------------------------------------------------
1028 subroutine gen_dist_constr
1029 C Generate CA distance constraints.
1030 implicit real*8 (a-h,o-z)
1031 include 'DIMENSIONS'
1032 include 'COMMON.IOUNITS'
1033 include 'COMMON.GEO'
1034 include 'COMMON.VAR'
1035 include 'COMMON.INTERACT'
1036 include 'COMMON.LOCAL'
1037 include 'COMMON.NAMES'
1038 include 'COMMON.CHAIN'
1039 include 'COMMON.FFIELD'
1040 include 'COMMON.SBRIDGE'
1041 include 'COMMON.HEADER'
1042 include 'COMMON.CONTROL'
1043 c include 'COMMON.DBASE'
1044 c include 'COMMON.THREAD'
1045 include 'COMMON.TIME1'
1046 dimension itype_pdb(maxres)
1047 common /pizda/ itype_pdb
1049 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1050 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1051 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1053 do i=nstart_sup,nstart_sup+nsup-1
1054 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1055 cd & ' seq_pdb', restyp(itype_pdb(i))
1056 do j=i+2,nstart_sup+nsup-1
1058 ihpb(nhpb)=i+nstart_seq-nstart_sup
1059 jhpb(nhpb)=j+nstart_seq-nstart_sup
1061 dhpb(nhpb)=dist(i,j)
1064 cd write (iout,'(a)') 'Distance constraints:'
1069 cd if (ii.gt.nres) then
1074 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1075 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1076 cd & dhpb(i),forcon(i)
1080 c----------------------------------------------------------------------------
1082 implicit real*8 (a-h,o-z)
1083 include 'DIMENSIONS'
1084 include 'COMMON.IOUNITS'
1085 include 'COMMON.GEO'
1086 include 'COMMON.CSA'
1087 include 'COMMON.BANK'
1088 include 'COMMON.CONTROL'
1090 character*620 mcmcard
1091 call card_concat(mcmcard)
1093 call readi(mcmcard,'NCONF',nconf,50)
1094 call readi(mcmcard,'NADD',nadd,0)
1095 call readi(mcmcard,'JSTART',jstart,1)
1096 call readi(mcmcard,'JEND',jend,1)
1097 call readi(mcmcard,'NSTMAX',nstmax,500000)
1098 call readi(mcmcard,'N0',n0,1)
1099 call readi(mcmcard,'N1',n1,6)
1100 call readi(mcmcard,'N2',n2,4)
1101 call readi(mcmcard,'N3',n3,0)
1102 call readi(mcmcard,'N4',n4,0)
1103 call readi(mcmcard,'N5',n5,0)
1104 call readi(mcmcard,'N6',n6,10)
1105 call readi(mcmcard,'N7',n7,0)
1106 call readi(mcmcard,'N8',n8,0)
1107 call readi(mcmcard,'N9',n9,0)
1108 call readi(mcmcard,'N14',n14,0)
1109 call readi(mcmcard,'N15',n15,0)
1110 call readi(mcmcard,'N16',n16,0)
1111 call readi(mcmcard,'N17',n17,0)
1112 call readi(mcmcard,'N18',n18,0)
1114 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1116 call readi(mcmcard,'NDIFF',ndiff,2)
1117 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1118 call readi(mcmcard,'IS1',is1,1)
1119 call readi(mcmcard,'IS2',is2,8)
1120 call readi(mcmcard,'NRAN0',nran0,4)
1121 call readi(mcmcard,'NRAN1',nran1,2)
1122 call readi(mcmcard,'IRR',irr,1)
1123 call readi(mcmcard,'NSEED',nseed,20)
1124 call readi(mcmcard,'NTOTAL',ntotal,10000)
1125 call reada(mcmcard,'CUT1',cut1,2.0d0)
1126 call reada(mcmcard,'CUT2',cut2,5.0d0)
1127 call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1128 call readi(mcmcard,'ICMAX',icmax,1)
1129 call readi(mcmcard,'NBANKM',nbankm,400)
1130 call readi(mcmcard,'IUCUT',iucut,2)
1131 call readi(mcmcard,'IRESTART',irestart,0)
1132 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1135 call reada(mcmcard,'DELE',dele,20.0d0)
1136 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1137 call readi(mcmcard,'IREF',iref,0)
1138 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1139 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1140 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1141 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1142 write (iout,*) "NCONF_IN",nconf_in
1143 tm_score=(index(mcmcard,'TMSCORE').gt.0)
1144 if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1145 & " for torsional angles"
1149 subroutine read_minim
1150 implicit real*8 (a-h,o-z)
1151 include 'DIMENSIONS'
1152 include 'COMMON.MINIM'
1153 include 'COMMON.IOUNITS'
1155 character*320 minimcard
1156 call card_concat(minimcard)
1157 call readi(minimcard,'MAXMIN',maxmin,2000)
1158 call readi(minimcard,'MAXFUN',maxfun,5000)
1159 call readi(minimcard,'MINMIN',minmin,maxmin)
1160 call readi(minimcard,'MINFUN',minfun,maxmin)
1161 call reada(minimcard,'TOLF',tolf,1.0D-2)
1162 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1163 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1164 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1165 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1166 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1167 & 'Options in energy minimization:'
1168 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1169 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1170 & 'MinMin:',MinMin,' MinFun:',MinFun,
1171 & ' TolF:',TolF,' RTolF:',RTolF
1174 c----------------------------------------------------------------------------
1175 subroutine read_angles(kanal,*)
1176 implicit real*8 (a-h,o-z)
1177 include 'DIMENSIONS'
1178 include 'COMMON.GEO'
1179 include 'COMMON.VAR'
1180 include 'COMMON.CHAIN'
1181 include 'COMMON.IOUNITS'
1182 include 'COMMON.CONTROL'
1183 c Read angles from input
1185 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1186 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1187 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1188 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1191 c 9/7/01 avoid 180 deg valence angle
1192 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1194 theta(i)=deg2rad*theta(i)
1195 phi(i)=deg2rad*phi(i)
1196 alph(i)=deg2rad*alph(i)
1197 omeg(i)=deg2rad*omeg(i)
1202 c----------------------------------------------------------------------------
1203 subroutine reada(rekord,lancuch,wartosc,default)
1205 character*(*) rekord,lancuch
1206 double precision wartosc,default
1209 iread=index(rekord,lancuch)
1210 if (iread.eq.0) then
1214 iread=iread+ilen(lancuch)+1
1215 read (rekord(iread:),*,err=10,end=10) wartosc
1220 c----------------------------------------------------------------------------
1221 subroutine readi(rekord,lancuch,wartosc,default)
1223 character*(*) rekord,lancuch
1224 integer wartosc,default
1227 iread=index(rekord,lancuch)
1228 if (iread.eq.0) then
1232 iread=iread+ilen(lancuch)+1
1233 read (rekord(iread:),*,err=10,end=10) wartosc
1238 c----------------------------------------------------------------------------
1239 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1242 integer tablica(dim),default
1243 character*(*) rekord,lancuch
1250 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1251 if (iread.eq.0) return
1252 iread=iread+ilen(lancuch)+1
1253 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1256 c----------------------------------------------------------------------------
1257 subroutine multreada(rekord,lancuch,tablica,dim,default)
1260 double precision tablica(dim),default
1261 character*(*) rekord,lancuch
1268 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1269 if (iread.eq.0) return
1270 iread=iread+ilen(lancuch)+1
1271 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1274 c----------------------------------------------------------------------------
1275 subroutine openunits
1276 implicit real*8 (a-h,o-z)
1277 include 'DIMENSIONS'
1280 character*16 form,nodename
1283 include 'COMMON.SETUP'
1284 include 'COMMON.IOUNITS'
1285 c include 'COMMON.MD'
1286 include 'COMMON.CONTROL'
1287 integer lenpre,lenpot,ilen,lentmp
1289 character*3 out1file_text,ucase
1292 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1293 call getenv_loc("PREFIX",prefix)
1295 call getenv_loc("POT",pot)
1296 call getenv_loc("DIRTMP",tmpdir)
1297 call getenv_loc("CURDIR",curdir)
1298 call getenv_loc("OUT1FILE",out1file_text)
1299 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1300 out1file_text=ucase(out1file_text)
1301 if (out1file_text(1:1).eq."Y") then
1304 out1file=fg_rank.gt.0
1309 if (lentmp.gt.0) then
1310 write (*,'(80(1h!))')
1311 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1312 write (*,'(80(1h!))')
1313 write (*,*)"All output files will be on node /tmp directory."
1315 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1316 if (me.eq.king) then
1317 write (*,*) "The master node is ",nodename
1318 else if (fg_rank.eq.0) then
1319 write (*,*) "I am the CG slave node ",nodename
1321 write (*,*) "I am the FG slave node ",nodename
1324 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1325 lenpre = lentmp+lenpre+1
1327 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1328 C Get the names and open the input files
1329 #if defined(WINIFL) || defined(WINPGI)
1330 open(1,file=pref_orig(:ilen(pref_orig))//
1331 & '.inp',status='old',readonly,shared)
1332 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1333 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1334 C Get parameter filenames and open the parameter files.
1335 call getenv_loc('BONDPAR',bondname)
1336 open (ibond,file=bondname,status='old',readonly,shared)
1337 call getenv_loc('THETPAR',thetname)
1338 open (ithep,file=thetname,status='old',readonly,shared)
1340 call getenv_loc('THETPARPDB',thetname_pdb)
1341 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1343 call getenv_loc('ROTPAR',rotname)
1344 open (irotam,file=rotname,status='old',readonly,shared)
1346 call getenv_loc('ROTPARPDB',rotname_pdb)
1347 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1349 call getenv_loc('TORPAR',torname)
1350 open (itorp,file=torname,status='old',readonly,shared)
1351 call getenv_loc('TORDPAR',tordname)
1352 open (itordp,file=tordname,status='old',readonly,shared)
1353 call getenv_loc('FOURIER',fouriername)
1354 open (ifourier,file=fouriername,status='old',readonly,shared)
1355 call getenv_loc('ELEPAR',elename)
1356 open (ielep,file=elename,status='old',readonly,shared)
1357 call getenv_loc('SIDEPAR',sidename)
1358 open (isidep,file=sidename,status='old',readonly,shared)
1359 #elif (defined CRAY) || (defined AIX)
1360 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1362 c print *,"Processor",myrank," opened file 1"
1363 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1364 c print *,"Processor",myrank," opened file 9"
1365 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1366 C Get parameter filenames and open the parameter files.
1367 call getenv_loc('BONDPAR',bondname)
1368 open (ibond,file=bondname,status='old',action='read')
1369 c print *,"Processor",myrank," opened file IBOND"
1370 call getenv_loc('THETPAR',thetname)
1371 open (ithep,file=thetname,status='old',action='read')
1372 c print *,"Processor",myrank," opened file ITHEP"
1374 call getenv_loc('THETPARPDB',thetname_pdb)
1375 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1377 call getenv_loc('ROTPAR',rotname)
1378 open (irotam,file=rotname,status='old',action='read')
1379 c print *,"Processor",myrank," opened file IROTAM"
1381 call getenv_loc('ROTPARPDB',rotname_pdb)
1382 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1384 call getenv_loc('TORPAR',torname)
1385 open (itorp,file=torname,status='old',action='read')
1386 c print *,"Processor",myrank," opened file ITORP"
1387 call getenv_loc('TORDPAR',tordname)
1388 open (itordp,file=tordname,status='old',action='read')
1389 c print *,"Processor",myrank," opened file ITORDP"
1390 call getenv_loc('SCCORPAR',sccorname)
1391 open (isccor,file=sccorname,status='old',action='read')
1392 c print *,"Processor",myrank," opened file ISCCOR"
1393 call getenv_loc('FOURIER',fouriername)
1394 open (ifourier,file=fouriername,status='old',action='read')
1395 c print *,"Processor",myrank," opened file IFOURIER"
1396 call getenv_loc('ELEPAR',elename)
1397 open (ielep,file=elename,status='old',action='read')
1398 c print *,"Processor",myrank," opened file IELEP"
1399 call getenv_loc('SIDEPAR',sidename)
1400 open (isidep,file=sidename,status='old',action='read')
1401 c print *,"Processor",myrank," opened file ISIDEP"
1402 c print *,"Processor",myrank," opened parameter files"
1404 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1405 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1406 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1407 C Get parameter filenames and open the parameter files.
1408 call getenv_loc('BONDPAR',bondname)
1409 open (ibond,file=bondname,status='old')
1410 call getenv_loc('THETPAR',thetname)
1411 open (ithep,file=thetname,status='old')
1413 call getenv_loc('THETPARPDB',thetname_pdb)
1414 open (ithep_pdb,file=thetname_pdb,status='old')
1416 call getenv_loc('ROTPAR',rotname)
1417 open (irotam,file=rotname,status='old')
1419 call getenv_loc('ROTPARPDB',rotname_pdb)
1420 open (irotam_pdb,file=rotname_pdb,status='old')
1422 call getenv_loc('TORPAR',torname)
1423 open (itorp,file=torname,status='old')
1424 call getenv_loc('TORDPAR',tordname)
1425 open (itordp,file=tordname,status='old')
1426 call getenv_loc('SCCORPAR',sccorname)
1427 open (isccor,file=sccorname,status='old')
1428 call getenv_loc('FOURIER',fouriername)
1429 open (ifourier,file=fouriername,status='old')
1430 call getenv_loc('ELEPAR',elename)
1431 open (ielep,file=elename,status='old')
1432 call getenv_loc('SIDEPAR',sidename)
1433 open (isidep,file=sidename,status='old')
1435 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1437 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1438 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1439 C Get parameter filenames and open the parameter files.
1440 call getenv_loc('BONDPAR',bondname)
1441 open (ibond,file=bondname,status='old',readonly)
1442 call getenv_loc('THETPAR',thetname)
1443 open (ithep,file=thetname,status='old',readonly)
1445 call getenv_loc('THETPARPDB',thetname_pdb)
1446 print *,"thetname_pdb ",thetname_pdb
1447 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1448 print *,ithep_pdb," opened"
1450 call getenv_loc('ROTPAR',rotname)
1451 open (irotam,file=rotname,status='old',readonly)
1453 call getenv_loc('ROTPARPDB',rotname_pdb)
1454 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1456 call getenv_loc('TORPAR',torname)
1457 open (itorp,file=torname,status='old',readonly)
1458 call getenv_loc('TORDPAR',tordname)
1459 open (itordp,file=tordname,status='old',readonly)
1460 call getenv_loc('SCCORPAR',sccorname)
1461 open (isccor,file=sccorname,status='old',readonly)
1462 call getenv_loc('FOURIER',fouriername)
1463 open (ifourier,file=fouriername,status='old',readonly)
1464 call getenv_loc('ELEPAR',elename)
1465 open (ielep,file=elename,status='old',readonly)
1466 call getenv_loc('SIDEPAR',sidename)
1467 open (isidep,file=sidename,status='old',readonly)
1471 C 8/9/01 In the newest version SCp interaction constants are read from a file
1472 C Use -DOLDSCP to use hard-coded constants instead.
1474 call getenv_loc('SCPPAR',scpname)
1475 #if defined(WINIFL) || defined(WINPGI)
1476 open (iscpp,file=scpname,status='old',readonly,shared)
1477 #elif (defined CRAY) || (defined AIX)
1478 open (iscpp,file=scpname,status='old',action='read')
1480 open (iscpp,file=scpname,status='old')
1482 open (iscpp,file=scpname,status='old',readonly)
1485 call getenv_loc('PATTERN',patname)
1486 #if defined(WINIFL) || defined(WINPGI)
1487 open (icbase,file=patname,status='old',readonly,shared)
1488 #elif (defined CRAY) || (defined AIX)
1489 open (icbase,file=patname,status='old',action='read')
1491 open (icbase,file=patname,status='old')
1493 open (icbase,file=patname,status='old',readonly)
1496 C Open output file only for CG processes
1497 c print *,"Processor",myrank," fg_rank",fg_rank
1498 if (fg_rank.eq.0) then
1500 if (nodes.eq.1) then
1503 npos = dlog10(dfloat(nodes-1))+1
1505 if (npos.lt.3) npos=3
1506 write (liczba,'(i1)') npos
1507 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1509 write (liczba,form) me
1510 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1511 & liczba(:ilen(liczba))
1512 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1514 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1516 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1517 & liczba(:ilen(liczba))//'.mol2'
1518 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1519 & liczba(:ilen(liczba))//'.stat'
1521 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1522 & //liczba(:ilen(liczba))//'.stat')
1523 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1526 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1527 c & liczba(:ilen(liczba))//'.const'
1532 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1533 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1534 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1535 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1536 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1538 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1540 rest2name=prefix(:ilen(prefix))//'.rst'
1542 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1545 #if defined(AIX) || defined(PGI)
1546 if (me.eq.king .or. .not. out1file)
1547 & open(iout,file=outname,status='unknown')
1549 if (fg_rank.gt.0) then
1550 write (liczba,'(i3.3)') myrank/nfgtasks
1551 write (ll,'(bz,i3.3)') fg_rank
1552 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1557 open(igeom,file=intname,status='unknown',position='append')
1558 open(ipdb,file=pdbname,status='unknown')
1559 open(imol2,file=mol2name,status='unknown')
1560 open(istat,file=statname,status='unknown',position='append')
1562 c1out open(iout,file=outname,status='unknown')
1565 if (me.eq.king .or. .not.out1file)
1566 & open(iout,file=outname,status='unknown')
1568 if (fg_rank.gt.0) then
1569 write (liczba,'(i3.3)') myrank/nfgtasks
1570 write (ll,'(bz,i3.3)') fg_rank
1571 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1576 open(igeom,file=intname,status='unknown',access='append')
1577 open(ipdb,file=pdbname,status='unknown')
1578 open(imol2,file=mol2name,status='unknown')
1579 open(istat,file=statname,status='unknown',access='append')
1581 c1out open(iout,file=outname,status='unknown')
1584 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1585 csa_seed=prefix(:lenpre)//'.CSA.seed'
1586 csa_history=prefix(:lenpre)//'.CSA.history'
1587 csa_bank=prefix(:lenpre)//'.CSA.bank'
1588 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1589 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1590 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1591 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1592 csa_int=prefix(:lenpre)//'.int'
1593 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1594 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1595 csa_in=prefix(:lenpre)//'.CSA.in'
1596 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1599 write (iout,'(80(1h-))')
1600 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1601 write (iout,'(80(1h-))')
1602 write (iout,*) "Input file : ",
1603 & pref_orig(:ilen(pref_orig))//'.inp'
1604 write (iout,*) "Output file : ",
1605 & outname(:ilen(outname))
1607 write (iout,*) "Sidechain potential file : ",
1608 & sidename(:ilen(sidename))
1610 write (iout,*) "SCp potential file : ",
1611 & scpname(:ilen(scpname))
1613 write (iout,*) "Electrostatic potential file : ",
1614 & elename(:ilen(elename))
1615 write (iout,*) "Cumulant coefficient file : ",
1616 & fouriername(:ilen(fouriername))
1617 write (iout,*) "Torsional parameter file : ",
1618 & torname(:ilen(torname))
1619 write (iout,*) "Double torsional parameter file : ",
1620 & tordname(:ilen(tordname))
1621 write (iout,*) "SCCOR parameter file : ",
1622 & sccorname(:ilen(sccorname))
1623 write (iout,*) "Bond & inertia constant file : ",
1624 & bondname(:ilen(bondname))
1625 write (iout,*) "Bending parameter file : ",
1626 & thetname(:ilen(thetname))
1627 write (iout,*) "Rotamer parameter file : ",
1628 & rotname(:ilen(rotname))
1629 write (iout,*) "Threading database : ",
1630 & patname(:ilen(patname))
1632 &write (iout,*)" DIRTMP : ",
1634 write (iout,'(80(1h-))')
1638 c----------------------------------------------------------------------------
1639 subroutine card_concat(card)
1640 implicit real*8 (a-h,o-z)
1641 include 'DIMENSIONS'
1642 include 'COMMON.IOUNITS'
1644 character*80 karta,ucase
1646 read (inp,'(a)') karta
1649 do while (karta(80:80).eq.'&')
1650 card=card(:ilen(card)+1)//karta(:79)
1651 read (inp,'(a)') karta
1654 card=card(:ilen(card)+1)//karta
1658 subroutine read_dist_constr
1659 implicit real*8 (a-h,o-z)
1660 include 'DIMENSIONS'
1664 include 'COMMON.SETUP'
1665 include 'COMMON.CONTROL'
1666 include 'COMMON.CHAIN'
1667 include 'COMMON.IOUNITS'
1668 include 'COMMON.SBRIDGE'
1669 integer ifrag_(2,100),ipair_(2,100)
1670 double precision wfrag_(100),wpair_(100)
1671 character*500 controlcard
1672 write (iout,*) "Calling read_dist_constr"
1673 write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1675 call card_concat(controlcard)
1676 call readi(controlcard,"NFRAG",nfrag_,0)
1677 call readi(controlcard,"NPAIR",npair_,0)
1678 call readi(controlcard,"NDIST",ndist_,0)
1679 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1680 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1681 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1682 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1683 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1684 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1685 c write (iout,*) "IFRAG"
1687 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1689 c write (iout,*) "IPAIR"
1691 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1695 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1696 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1697 & ifrag_(2,i)=nstart_sup+nsup-1
1698 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1700 if (wfrag_(i).gt.0.0d0) then
1701 do j=ifrag_(1,i),ifrag_(2,i)-1
1702 do k=j+1,ifrag_(2,i)
1703 write (iout,*) "j",j," k",k
1705 if (constr_dist.eq.1) then
1710 forcon(nhpb)=wfrag_(i)
1711 else if (constr_dist.eq.2) then
1712 if (ddjk.le.dist_cut) then
1717 forcon(nhpb)=wfrag_(i)
1724 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1727 if (.not.out1file .or. me.eq.king)
1728 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1729 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1731 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1732 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1739 if (wpair_(i).gt.0.0d0) then
1747 do j=ifrag_(1,ii),ifrag_(2,ii)
1748 do k=ifrag_(1,jj),ifrag_(2,jj)
1752 forcon(nhpb)=wpair_(i)
1753 dhpb(nhpb)=dist(j,k)
1755 if (.not.out1file .or. me.eq.king)
1756 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1757 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1759 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1760 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1767 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1768 if (forcon(nhpb+1).gt.0.0d0) then
1770 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1772 if (.not.out1file .or. me.eq.king)
1773 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1774 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1776 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1777 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1784 c-------------------------------------------------------------------------------
1786 subroutine flush(iu)
1791 subroutine flush(iu)
1796 c------------------------------------------------------------------------------
1797 subroutine copy_to_tmp(source)
1798 include "DIMENSIONS"
1799 include "COMMON.IOUNITS"
1800 character*(*) source
1801 character* 256 tmpfile
1805 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1806 inquire(file=tmpfile,exist=ex)
1808 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1809 & " to temporary directory..."
1810 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1811 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1815 c------------------------------------------------------------------------------
1816 subroutine move_from_tmp(source)
1817 include "DIMENSIONS"
1818 include "COMMON.IOUNITS"
1819 character*(*) source
1822 write (*,*) "Moving ",source(:ilen(source)),
1823 & " from temporary directory to working directory"
1824 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1825 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1828 c------------------------------------------------------------------------------
1829 subroutine random_init(seed)
1831 C Initialize random number generator
1833 implicit real*8 (a-h,o-z)
1834 include 'DIMENSIONS'
1840 logical OKRandom, prng_restart
1842 integer iseed_array(4)
1844 include 'COMMON.IOUNITS'
1845 include 'COMMON.TIME1'
1846 c include 'COMMON.THREAD'
1847 include 'COMMON.SBRIDGE'
1848 include 'COMMON.CONTROL'
1849 include 'COMMON.MCM'
1850 c include 'COMMON.MAP'
1851 include 'COMMON.HEADER'
1852 c include 'COMMON.CSA'
1853 include 'COMMON.CHAIN'
1854 c include 'COMMON.MUCA'
1855 c include 'COMMON.MD'
1856 include 'COMMON.FFIELD'
1857 include 'COMMON.SETUP'
1858 iseed=-dint(dabs(seed))
1859 if (iseed.eq.0) then
1860 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1861 & 'Random seed undefined. The program will stop.'
1862 write (*,'(/80(1h*)/20x,a/80(1h*))')
1863 & 'Random seed undefined. The program will stop.'
1865 call mpi_finalize(mpi_comm_world,ierr)
1867 stop 'Bad random seed.'
1870 if (fg_rank.eq.0) then
1874 if(me.eq.king .or. .not. out1file)
1875 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1876 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1877 OKRandom = prng_restart(me,iseedi8)
1880 tmp=65536.0d0**(4-i)
1881 iseed_array(i) = dint(seed/tmp)
1882 seed=seed-iseed_array(i)*tmp
1884 if(me.eq.king .or. .not. out1file)
1885 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1886 & (iseed_array(i),i=1,4)
1887 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1888 & (iseed_array(i),i=1,4)
1889 OKRandom = prng_restart(me,iseed_array)
1892 r1=ran_number(0.0D0,1.0D0)
1893 if(me.eq.king .or. .not. out1file)
1894 & write (iout,*) 'ran_num',r1
1895 if (r1.lt.0.0d0) OKRandom=.false.
1897 if (.not.OKRandom) then
1898 write (iout,*) 'PRNG IS NOT WORKING!!!'
1899 print *,'PRNG IS NOT WORKING!!!'
1902 call mpi_abort(mpi_comm_world,error_msg,ierr)
1905 write (iout,*) 'too many processors for parallel prng'
1906 write (*,*) 'too many processors for parallel prng'
1914 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)