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 C Following 2 lines for diagnostics; comment out if not needed
459 write (iout,*) "Before sideadd"
461 if(me.eq.king.or..not.out1file)
462 & write(iout,*)'Adding sidechains'
469 do while (fail.and.nsi.le.maxsi)
470 c call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
473 if(fail) write(iout,*)'Adding sidechain failed for res ',
474 & i,' after ',nsi,' trials'
478 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
480 C Following 2 lines for diagnostics; comment out if not needed
481 write (iout,*) "After sideadd"
485 if (indpdb.eq.0) then
486 C Read sequence if not taken from the pdb file.
488 c print *,'nres=',nres
489 if (iscode.gt.0) then
490 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
492 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
494 C Convert sequence to numeric code
496 itype(i)=rescode(i,sequence(i),iscode)
498 C Assign initial virtual bond lengths
504 vbld(i+nres)=dsc(iabs(itype(i)))
505 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
506 c write (iout,*) "i",i," itype",itype(i),
507 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
511 c print '(20i4)',(itype(i),i=1,nres)
514 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
516 if (itype(i).eq.21) then
520 else if (iabs(itype(i+1)).ne.20) then
522 else if (iabs(itype(i)).ne.20) then
529 if(me.eq.king.or..not.out1file)then
530 write (iout,*) "ITEL"
532 write (iout,*) i,itype(i),itel(i)
534 print *,'Call Read_Bridge.'
537 C 8/13/98 Set limits to generating the dihedral angles
542 read (inp,*) ndih_constr
543 if (ndih_constr.gt.0) then
545 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
546 if(me.eq.king.or..not.out1file)then
548 & 'There are',ndih_constr,' constraints on phi angles.'
550 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
554 phi0(i)=deg2rad*phi0(i)
555 drange(i)=deg2rad*drange(i)
557 if(me.eq.king.or..not.out1file)
558 & write (iout,*) 'FTORS',ftors
561 phibound(1,ii) = phi0(i)-drange(i)
562 phibound(2,ii) = phi0(i)+drange(i)
569 write (iout,'(a)') 'Boundaries in phi angle sampling:'
571 write (iout,'(a3,i5,2f10.1)')
572 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
578 cd print *,'NNT=',NNT,' NCT=',NCT
579 if (itype(1).eq.21) nnt=2
580 if (itype(nres).eq.21) nct=nct-1
582 C Juyong:READ init_vars
583 C Initialize variables!
584 C Juyong:READ read_info
585 C READ fragment information!!
586 C both routines should be in dfa.F file!!
588 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
589 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
591 print*, 'init_dfa_vars finished!'
593 print*, 'read_dfa_info finished!'
600 if(me.eq.king.or..not.out1file)
601 & write (iout,'(a,i3)') 'nsup=',nsup
603 if (nsup.le.(nct-nnt+1)) then
604 do i=0,nct-nnt+1-nsup
605 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
611 & 'Error - sequences to be superposed do not match.'
614 do i=0,nsup-(nct-nnt+1)
615 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
617 nstart_sup=nstart_sup+i
623 & 'Error - sequences to be superposed do not match.'
626 if (nsup.eq.0) nsup=nct-nnt
627 if (nstart_sup.eq.0) nstart_sup=nnt
628 if (nstart_seq.eq.0) nstart_seq=nnt
629 if(me.eq.king.or..not.out1file)
630 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
631 & ' nstart_seq=',nstart_seq
633 c--- Zscore rms -------
634 if (nz_start.eq.0) nz_start=nnt
635 if (nz_end.eq.0 .and. nsup.gt.0) then
637 else if (nz_end.eq.0) then
640 if(me.eq.king.or..not.out1file)then
641 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
642 write (iout,*) 'IZ_SC=',iz_sc
644 c----------------------
647 if (.not.pdbref) then
648 call read_angles(inp,*38)
650 38 write (iout,'(a)') 'Error reading reference structure.'
652 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
653 stop 'Error reading reference structure'
657 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
666 call contact(.true.,ncont_ref,icont_ref,co)
668 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
670 if (constr_dist.gt.0) call read_dist_constr
671 c write (iout,*) "After read_dist_constr nhpb",nhpb
673 if(me.eq.king.or..not.out1file)
674 & write (iout,*) 'Contact order:',co
676 if(me.eq.king.or..not.out1file)
677 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
680 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
682 if(me.eq.king.or..not.out1file)
683 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
684 & icont_ref(1,i),' ',
685 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
689 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
690 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
691 & modecalc.ne.10) then
692 C If input structure hasn't been supplied from the PDB file read or generate
694 if (iranconf.eq.0 .and. .not. extconf) then
695 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
696 & write (iout,'(a)') 'Initial geometry will be read in.'
698 read(inp,'(8f10.5)',end=36,err=36)
699 & ((c(l,k),l=1,3),k=1,nres),
700 & ((c(l,k+nres),l=1,3),k=nnt,nct)
701 call int_from_cart1(.false.)
704 dc(j,i)=c(j,i+1)-c(j,i)
705 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
709 if (itype(i).ne.10) then
711 dc(j,i+nres)=c(j,i+nres)-c(j,i)
712 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
718 call read_angles(inp,*36)
721 36 write (iout,'(a)') 'Error reading angle file.'
723 call mpi_finalize( MPI_COMM_WORLD,IERR )
725 stop 'Error reading angle file.'
727 else if (extconf) then
728 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
729 & write (iout,'(a)') 'Extended chain initial geometry.'
731 theta(i)=90d0*deg2rad
737 alph(i)=110d0*deg2rad
740 omeg(i)=-120d0*deg2rad
741 if (itype(i).le.0) omeg(i)=-omeg(i)
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.
812 if (nthread.gt.0) then
816 if (me.eq.king .or. .not. out1file)
818 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
819 write (iout,'(/a,i3,a)')
820 & 'The chain contains',ns,' disulfide-bridging cysteines.'
821 write (iout,'(20i4)') (iss(i),i=1,ns)
822 write (iout,'(/a/)') 'Pre-formed links are:'
828 if (me.eq.king.or..not.out1file)
829 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
830 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
835 c if (i2ndstr.gt.0) call secstrp2dihc
836 c call geom_to_var(nvar,x)
837 c call etotal(energia(0))
838 c call enerprint(energia(0))
839 c call briefout(0,etot)
841 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
842 cd write (iout,'(a)') 'Variable list:'
843 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
845 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
846 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
847 & 'Processor',myrank,': end reading molecular data.'
851 c--------------------------------------------------------------------------
852 logical function seq_comp(itypea,itypeb,length)
854 integer length,itypea(length),itypeb(length)
857 if (itypea(i).ne.itypeb(i)) then
865 c-----------------------------------------------------------------------------
866 subroutine read_bridge
867 C Read information about disulfide bridges.
868 implicit real*8 (a-h,o-z)
873 include 'COMMON.IOUNITS'
876 include 'COMMON.INTERACT'
877 include 'COMMON.LOCAL'
878 include 'COMMON.NAMES'
879 include 'COMMON.CHAIN'
880 include 'COMMON.FFIELD'
881 include 'COMMON.SBRIDGE'
882 include 'COMMON.HEADER'
883 include 'COMMON.CONTROL'
884 c include 'COMMON.DBASE'
885 c include 'COMMON.THREAD'
886 include 'COMMON.TIME1'
887 include 'COMMON.SETUP'
888 C Read bridging residues.
889 read (inp,*) ns,(iss(i),i=1,ns)
891 if(me.eq.king.or..not.out1file)
892 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
893 C Check whether the specified bridging residues are cystines.
895 if (itype(iss(i)).ne.1) then
896 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
897 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
898 & ' can form a disulfide bridge?!!!'
899 write (*,'(2a,i3,a)')
900 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
901 & ' can form a disulfide bridge?!!!'
903 call MPI_Finalize(MPI_COMM_WORLD,ierror)
908 C Read preformed bridges.
910 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
911 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
914 C Check if the residues involved in bridges are in the specified list of
918 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
919 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
920 write (iout,'(a,i3,a)') 'Disulfide pair',i,
921 & ' contains residues present in other pairs.'
922 write (*,'(a,i3,a)') 'Disulfide pair',i,
923 & ' contains residues present in other pairs.'
925 call MPI_Finalize(MPI_COMM_WORLD,ierror)
931 if (ihpb(i).eq.iss(j)) goto 10
933 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
936 if (jhpb(i).eq.iss(j)) goto 20
938 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
951 c----------------------------------------------------------------------------
952 subroutine read_x(kanal,*)
953 implicit real*8 (a-h,o-z)
957 include 'COMMON.CHAIN'
958 include 'COMMON.IOUNITS'
959 include 'COMMON.CONTROL'
960 include 'COMMON.LOCAL'
961 include 'COMMON.INTERACT'
962 c Read coordinates from input
964 read(kanal,'(8f10.5)',end=10,err=10)
965 & ((c(l,k),l=1,3),k=1,nres),
966 & ((c(l,k+nres),l=1,3),k=nnt,nct)
969 c(j,2*nres)=c(j,nres)
971 call int_from_cart1(.false.)
974 dc(j,i)=c(j,i+1)-c(j,i)
975 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
979 if (itype(i).ne.10) then
981 dc(j,i+nres)=c(j,i+nres)-c(j,i)
982 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
990 c------------------------------------------------------------------------------
992 implicit real*8 (a-h,o-z)
994 include 'COMMON.IOUNITS'
997 include 'COMMON.INTERACT'
998 include 'COMMON.LOCAL'
999 include 'COMMON.NAMES'
1000 include 'COMMON.CHAIN'
1001 include 'COMMON.FFIELD'
1002 include 'COMMON.SBRIDGE'
1003 include 'COMMON.HEADER'
1004 include 'COMMON.CONTROL'
1005 c include 'COMMON.DBASE'
1006 c include 'COMMON.THREAD'
1007 include 'COMMON.TIME1'
1008 C Set up variable list.
1014 if (itype(i).ne.10) then
1016 ialph(i,1)=nvar+nside
1020 if (indphi.gt.0) then
1022 else if (indback.gt.0) then
1027 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1030 c----------------------------------------------------------------------------
1031 subroutine gen_dist_constr
1032 C Generate CA distance constraints.
1033 implicit real*8 (a-h,o-z)
1034 include 'DIMENSIONS'
1035 include 'COMMON.IOUNITS'
1036 include 'COMMON.GEO'
1037 include 'COMMON.VAR'
1038 include 'COMMON.INTERACT'
1039 include 'COMMON.LOCAL'
1040 include 'COMMON.NAMES'
1041 include 'COMMON.CHAIN'
1042 include 'COMMON.FFIELD'
1043 include 'COMMON.SBRIDGE'
1044 include 'COMMON.HEADER'
1045 include 'COMMON.CONTROL'
1046 c include 'COMMON.DBASE'
1047 c include 'COMMON.THREAD'
1048 include 'COMMON.TIME1'
1049 dimension itype_pdb(maxres)
1050 common /pizda/ itype_pdb
1052 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1053 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1054 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1056 do i=nstart_sup,nstart_sup+nsup-1
1057 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1058 cd & ' seq_pdb', restyp(itype_pdb(i))
1059 do j=i+2,nstart_sup+nsup-1
1061 ihpb(nhpb)=i+nstart_seq-nstart_sup
1062 jhpb(nhpb)=j+nstart_seq-nstart_sup
1064 dhpb(nhpb)=dist(i,j)
1067 cd write (iout,'(a)') 'Distance constraints:'
1072 cd if (ii.gt.nres) then
1077 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1078 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1079 cd & dhpb(i),forcon(i)
1083 c----------------------------------------------------------------------------
1085 implicit real*8 (a-h,o-z)
1086 include 'DIMENSIONS'
1087 include 'COMMON.IOUNITS'
1088 include 'COMMON.GEO'
1089 include 'COMMON.CSA'
1090 include 'COMMON.BANK'
1091 include 'COMMON.CONTROL'
1093 character*620 mcmcard
1094 call card_concat(mcmcard)
1096 call readi(mcmcard,'NCONF',nconf,50)
1097 call readi(mcmcard,'NADD',nadd,0)
1098 call readi(mcmcard,'JSTART',jstart,1)
1099 call readi(mcmcard,'JEND',jend,1)
1100 call readi(mcmcard,'NSTMAX',nstmax,500000)
1101 call readi(mcmcard,'N0',n0,1)
1102 call readi(mcmcard,'N1',n1,6)
1103 call readi(mcmcard,'N2',n2,4)
1104 call readi(mcmcard,'N3',n3,0)
1105 call readi(mcmcard,'N4',n4,0)
1106 call readi(mcmcard,'N5',n5,0)
1107 call readi(mcmcard,'N6',n6,10)
1108 call readi(mcmcard,'N7',n7,0)
1109 call readi(mcmcard,'N8',n8,0)
1110 call readi(mcmcard,'N9',n9,0)
1111 call readi(mcmcard,'N14',n14,0)
1112 call readi(mcmcard,'N15',n15,0)
1113 call readi(mcmcard,'N16',n16,0)
1114 call readi(mcmcard,'N17',n17,0)
1115 call readi(mcmcard,'N18',n18,0)
1117 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1119 call readi(mcmcard,'NDIFF',ndiff,2)
1120 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1121 call readi(mcmcard,'IS1',is1,1)
1122 call readi(mcmcard,'IS2',is2,8)
1123 call readi(mcmcard,'NRAN0',nran0,4)
1124 call readi(mcmcard,'NRAN1',nran1,2)
1125 call readi(mcmcard,'IRR',irr,1)
1126 call readi(mcmcard,'NSEED',nseed,20)
1127 call readi(mcmcard,'NTOTAL',ntotal,10000)
1128 call reada(mcmcard,'CUT1',cut1,2.0d0)
1129 call reada(mcmcard,'CUT2',cut2,5.0d0)
1130 call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1131 call readi(mcmcard,'ICMAX',icmax,1)
1132 call readi(mcmcard,'NBANKM',nbankm,400)
1133 call readi(mcmcard,'IUCUT',iucut,2)
1134 call readi(mcmcard,'IRESTART',irestart,0)
1135 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1138 call reada(mcmcard,'DELE',dele,20.0d0)
1139 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1140 call readi(mcmcard,'IREF',iref,0)
1141 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1142 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1143 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1144 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1145 write (iout,*) "NCONF_IN",nconf_in
1146 tm_score=(index(mcmcard,'TMSCORE').gt.0)
1147 if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1148 & " for torsional angles"
1152 subroutine read_minim
1153 implicit real*8 (a-h,o-z)
1154 include 'DIMENSIONS'
1155 include 'COMMON.MINIM'
1156 include 'COMMON.IOUNITS'
1158 character*320 minimcard
1159 call card_concat(minimcard)
1160 call readi(minimcard,'MAXMIN',maxmin,2000)
1161 call readi(minimcard,'MAXFUN',maxfun,5000)
1162 call readi(minimcard,'MINMIN',minmin,maxmin)
1163 call readi(minimcard,'MINFUN',minfun,maxmin)
1164 call reada(minimcard,'TOLF',tolf,1.0D-2)
1165 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1166 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1167 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1168 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1169 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1170 & 'Options in energy minimization:'
1171 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1172 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1173 & 'MinMin:',MinMin,' MinFun:',MinFun,
1174 & ' TolF:',TolF,' RTolF:',RTolF
1177 c----------------------------------------------------------------------------
1178 subroutine read_angles(kanal,*)
1179 implicit real*8 (a-h,o-z)
1180 include 'DIMENSIONS'
1181 include 'COMMON.GEO'
1182 include 'COMMON.VAR'
1183 include 'COMMON.CHAIN'
1184 include 'COMMON.IOUNITS'
1185 include 'COMMON.CONTROL'
1186 c Read angles from input
1188 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1189 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1190 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1191 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1194 c 9/7/01 avoid 180 deg valence angle
1195 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1197 theta(i)=deg2rad*theta(i)
1198 phi(i)=deg2rad*phi(i)
1199 alph(i)=deg2rad*alph(i)
1200 omeg(i)=deg2rad*omeg(i)
1205 c----------------------------------------------------------------------------
1206 subroutine reada(rekord,lancuch,wartosc,default)
1208 character*(*) rekord,lancuch
1209 double precision wartosc,default
1212 iread=index(rekord,lancuch)
1213 if (iread.eq.0) then
1217 iread=iread+ilen(lancuch)+1
1218 read (rekord(iread:),*,err=10,end=10) wartosc
1223 c----------------------------------------------------------------------------
1224 subroutine readi(rekord,lancuch,wartosc,default)
1226 character*(*) rekord,lancuch
1227 integer wartosc,default
1230 iread=index(rekord,lancuch)
1231 if (iread.eq.0) then
1235 iread=iread+ilen(lancuch)+1
1236 read (rekord(iread:),*,err=10,end=10) wartosc
1241 c----------------------------------------------------------------------------
1242 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1245 integer tablica(dim),default
1246 character*(*) rekord,lancuch
1253 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1254 if (iread.eq.0) return
1255 iread=iread+ilen(lancuch)+1
1256 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1259 c----------------------------------------------------------------------------
1260 subroutine multreada(rekord,lancuch,tablica,dim,default)
1263 double precision tablica(dim),default
1264 character*(*) rekord,lancuch
1271 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1272 if (iread.eq.0) return
1273 iread=iread+ilen(lancuch)+1
1274 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1277 c----------------------------------------------------------------------------
1278 subroutine openunits
1279 implicit real*8 (a-h,o-z)
1280 include 'DIMENSIONS'
1283 character*16 form,nodename
1286 include 'COMMON.SETUP'
1287 include 'COMMON.IOUNITS'
1288 c include 'COMMON.MD'
1289 include 'COMMON.CONTROL'
1290 integer lenpre,lenpot,ilen,lentmp
1292 character*3 out1file_text,ucase
1295 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1296 call getenv_loc("PREFIX",prefix)
1298 call getenv_loc("POT",pot)
1299 call getenv_loc("DIRTMP",tmpdir)
1300 call getenv_loc("CURDIR",curdir)
1301 call getenv_loc("OUT1FILE",out1file_text)
1302 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1303 out1file_text=ucase(out1file_text)
1304 if (out1file_text(1:1).eq."Y") then
1307 out1file=fg_rank.gt.0
1312 if (lentmp.gt.0) then
1313 write (*,'(80(1h!))')
1314 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1315 write (*,'(80(1h!))')
1316 write (*,*)"All output files will be on node /tmp directory."
1318 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1319 if (me.eq.king) then
1320 write (*,*) "The master node is ",nodename
1321 else if (fg_rank.eq.0) then
1322 write (*,*) "I am the CG slave node ",nodename
1324 write (*,*) "I am the FG slave node ",nodename
1327 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1328 lenpre = lentmp+lenpre+1
1330 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1331 C Get the names and open the input files
1332 #if defined(WINIFL) || defined(WINPGI)
1333 open(1,file=pref_orig(:ilen(pref_orig))//
1334 & '.inp',status='old',readonly,shared)
1335 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1336 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1337 C Get parameter filenames and open the parameter files.
1338 call getenv_loc('BONDPAR',bondname)
1339 open (ibond,file=bondname,status='old',readonly,shared)
1340 call getenv_loc('THETPAR',thetname)
1341 open (ithep,file=thetname,status='old',readonly,shared)
1343 call getenv_loc('THETPARPDB',thetname_pdb)
1344 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1346 call getenv_loc('ROTPAR',rotname)
1347 open (irotam,file=rotname,status='old',readonly,shared)
1349 call getenv_loc('ROTPARPDB',rotname_pdb)
1350 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1352 call getenv_loc('TORPAR',torname)
1353 open (itorp,file=torname,status='old',readonly,shared)
1354 call getenv_loc('TORDPAR',tordname)
1355 open (itordp,file=tordname,status='old',readonly,shared)
1356 call getenv_loc('FOURIER',fouriername)
1357 open (ifourier,file=fouriername,status='old',readonly,shared)
1358 call getenv_loc('ELEPAR',elename)
1359 open (ielep,file=elename,status='old',readonly,shared)
1360 call getenv_loc('SIDEPAR',sidename)
1361 open (isidep,file=sidename,status='old',readonly,shared)
1362 #elif (defined CRAY) || (defined AIX)
1363 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1365 c print *,"Processor",myrank," opened file 1"
1366 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1367 c print *,"Processor",myrank," opened file 9"
1368 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1369 C Get parameter filenames and open the parameter files.
1370 call getenv_loc('BONDPAR',bondname)
1371 open (ibond,file=bondname,status='old',action='read')
1372 c print *,"Processor",myrank," opened file IBOND"
1373 call getenv_loc('THETPAR',thetname)
1374 open (ithep,file=thetname,status='old',action='read')
1375 c print *,"Processor",myrank," opened file ITHEP"
1377 call getenv_loc('THETPARPDB',thetname_pdb)
1378 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1380 call getenv_loc('ROTPAR',rotname)
1381 open (irotam,file=rotname,status='old',action='read')
1382 c print *,"Processor",myrank," opened file IROTAM"
1384 call getenv_loc('ROTPARPDB',rotname_pdb)
1385 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1387 call getenv_loc('TORPAR',torname)
1388 open (itorp,file=torname,status='old',action='read')
1389 c print *,"Processor",myrank," opened file ITORP"
1390 call getenv_loc('TORDPAR',tordname)
1391 open (itordp,file=tordname,status='old',action='read')
1392 c print *,"Processor",myrank," opened file ITORDP"
1393 call getenv_loc('SCCORPAR',sccorname)
1394 open (isccor,file=sccorname,status='old',action='read')
1395 c print *,"Processor",myrank," opened file ISCCOR"
1396 call getenv_loc('FOURIER',fouriername)
1397 open (ifourier,file=fouriername,status='old',action='read')
1398 c print *,"Processor",myrank," opened file IFOURIER"
1399 call getenv_loc('ELEPAR',elename)
1400 open (ielep,file=elename,status='old',action='read')
1401 c print *,"Processor",myrank," opened file IELEP"
1402 call getenv_loc('SIDEPAR',sidename)
1403 open (isidep,file=sidename,status='old',action='read')
1404 c print *,"Processor",myrank," opened file ISIDEP"
1405 c print *,"Processor",myrank," opened parameter files"
1407 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1408 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1409 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1410 C Get parameter filenames and open the parameter files.
1411 call getenv_loc('BONDPAR',bondname)
1412 open (ibond,file=bondname,status='old')
1413 call getenv_loc('THETPAR',thetname)
1414 open (ithep,file=thetname,status='old')
1416 call getenv_loc('THETPARPDB',thetname_pdb)
1417 open (ithep_pdb,file=thetname_pdb,status='old')
1419 call getenv_loc('ROTPAR',rotname)
1420 open (irotam,file=rotname,status='old')
1422 call getenv_loc('ROTPARPDB',rotname_pdb)
1423 open (irotam_pdb,file=rotname_pdb,status='old')
1425 call getenv_loc('TORPAR',torname)
1426 open (itorp,file=torname,status='old')
1427 call getenv_loc('TORDPAR',tordname)
1428 open (itordp,file=tordname,status='old')
1429 call getenv_loc('SCCORPAR',sccorname)
1430 open (isccor,file=sccorname,status='old')
1431 call getenv_loc('FOURIER',fouriername)
1432 open (ifourier,file=fouriername,status='old')
1433 call getenv_loc('ELEPAR',elename)
1434 open (ielep,file=elename,status='old')
1435 call getenv_loc('SIDEPAR',sidename)
1436 open (isidep,file=sidename,status='old')
1438 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1440 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1441 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1442 C Get parameter filenames and open the parameter files.
1443 call getenv_loc('BONDPAR',bondname)
1444 open (ibond,file=bondname,status='old',readonly)
1445 call getenv_loc('THETPAR',thetname)
1446 open (ithep,file=thetname,status='old',readonly)
1448 call getenv_loc('THETPARPDB',thetname_pdb)
1449 print *,"thetname_pdb ",thetname_pdb
1450 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1451 print *,ithep_pdb," opened"
1453 call getenv_loc('ROTPAR',rotname)
1454 open (irotam,file=rotname,status='old',readonly)
1456 call getenv_loc('ROTPARPDB',rotname_pdb)
1457 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1459 call getenv_loc('TORPAR',torname)
1460 open (itorp,file=torname,status='old',readonly)
1461 call getenv_loc('TORDPAR',tordname)
1462 open (itordp,file=tordname,status='old',readonly)
1463 call getenv_loc('SCCORPAR',sccorname)
1464 open (isccor,file=sccorname,status='old',readonly)
1465 call getenv_loc('FOURIER',fouriername)
1466 open (ifourier,file=fouriername,status='old',readonly)
1467 call getenv_loc('ELEPAR',elename)
1468 open (ielep,file=elename,status='old',readonly)
1469 call getenv_loc('SIDEPAR',sidename)
1470 open (isidep,file=sidename,status='old',readonly)
1474 C 8/9/01 In the newest version SCp interaction constants are read from a file
1475 C Use -DOLDSCP to use hard-coded constants instead.
1477 call getenv_loc('SCPPAR',scpname)
1478 #if defined(WINIFL) || defined(WINPGI)
1479 open (iscpp,file=scpname,status='old',readonly,shared)
1480 #elif (defined CRAY) || (defined AIX)
1481 open (iscpp,file=scpname,status='old',action='read')
1483 open (iscpp,file=scpname,status='old')
1485 open (iscpp,file=scpname,status='old',readonly)
1488 call getenv_loc('PATTERN',patname)
1489 #if defined(WINIFL) || defined(WINPGI)
1490 open (icbase,file=patname,status='old',readonly,shared)
1491 #elif (defined CRAY) || (defined AIX)
1492 open (icbase,file=patname,status='old',action='read')
1494 open (icbase,file=patname,status='old')
1496 open (icbase,file=patname,status='old',readonly)
1499 C Open output file only for CG processes
1500 c print *,"Processor",myrank," fg_rank",fg_rank
1501 if (fg_rank.eq.0) then
1503 if (nodes.eq.1) then
1506 npos = dlog10(dfloat(nodes-1))+1
1508 if (npos.lt.3) npos=3
1509 write (liczba,'(i1)') npos
1510 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1512 write (liczba,form) me
1513 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1514 & liczba(:ilen(liczba))
1515 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1517 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1519 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1520 & liczba(:ilen(liczba))//'.mol2'
1521 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1522 & liczba(:ilen(liczba))//'.stat'
1524 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1525 & //liczba(:ilen(liczba))//'.stat')
1526 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1529 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1530 c & liczba(:ilen(liczba))//'.const'
1535 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1536 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1537 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1538 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1539 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1541 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1543 rest2name=prefix(:ilen(prefix))//'.rst'
1545 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1548 #if defined(AIX) || defined(PGI)
1549 if (me.eq.king .or. .not. out1file)
1550 & open(iout,file=outname,status='unknown')
1552 if (fg_rank.gt.0) then
1553 write (liczba,'(i3.3)') myrank/nfgtasks
1554 write (ll,'(bz,i3.3)') fg_rank
1555 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1560 open(igeom,file=intname,status='unknown',position='append')
1561 open(ipdb,file=pdbname,status='unknown')
1562 open(imol2,file=mol2name,status='unknown')
1563 open(istat,file=statname,status='unknown',position='append')
1565 c1out open(iout,file=outname,status='unknown')
1568 if (me.eq.king .or. .not.out1file)
1569 & open(iout,file=outname,status='unknown')
1571 if (fg_rank.gt.0) then
1572 write (liczba,'(i3.3)') myrank/nfgtasks
1573 write (ll,'(bz,i3.3)') fg_rank
1574 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1579 open(igeom,file=intname,status='unknown',access='append')
1580 open(ipdb,file=pdbname,status='unknown')
1581 open(imol2,file=mol2name,status='unknown')
1582 open(istat,file=statname,status='unknown',access='append')
1584 c1out open(iout,file=outname,status='unknown')
1587 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1588 csa_seed=prefix(:lenpre)//'.CSA.seed'
1589 csa_history=prefix(:lenpre)//'.CSA.history'
1590 csa_bank=prefix(:lenpre)//'.CSA.bank'
1591 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1592 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1593 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1594 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1595 csa_int=prefix(:lenpre)//'.int'
1596 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1597 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1598 csa_in=prefix(:lenpre)//'.CSA.in'
1599 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1602 write (iout,'(80(1h-))')
1603 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1604 write (iout,'(80(1h-))')
1605 write (iout,*) "Input file : ",
1606 & pref_orig(:ilen(pref_orig))//'.inp'
1607 write (iout,*) "Output file : ",
1608 & outname(:ilen(outname))
1610 write (iout,*) "Sidechain potential file : ",
1611 & sidename(:ilen(sidename))
1613 write (iout,*) "SCp potential file : ",
1614 & scpname(:ilen(scpname))
1616 write (iout,*) "Electrostatic potential file : ",
1617 & elename(:ilen(elename))
1618 write (iout,*) "Cumulant coefficient file : ",
1619 & fouriername(:ilen(fouriername))
1620 write (iout,*) "Torsional parameter file : ",
1621 & torname(:ilen(torname))
1622 write (iout,*) "Double torsional parameter file : ",
1623 & tordname(:ilen(tordname))
1624 write (iout,*) "SCCOR parameter file : ",
1625 & sccorname(:ilen(sccorname))
1626 write (iout,*) "Bond & inertia constant file : ",
1627 & bondname(:ilen(bondname))
1628 write (iout,*) "Bending parameter file : ",
1629 & thetname(:ilen(thetname))
1630 write (iout,*) "Rotamer parameter file : ",
1631 & rotname(:ilen(rotname))
1632 write (iout,*) "Threading database : ",
1633 & patname(:ilen(patname))
1635 &write (iout,*)" DIRTMP : ",
1637 write (iout,'(80(1h-))')
1641 c----------------------------------------------------------------------------
1642 subroutine card_concat(card)
1643 implicit real*8 (a-h,o-z)
1644 include 'DIMENSIONS'
1645 include 'COMMON.IOUNITS'
1647 character*80 karta,ucase
1649 read (inp,'(a)') karta
1652 do while (karta(80:80).eq.'&')
1653 card=card(:ilen(card)+1)//karta(:79)
1654 read (inp,'(a)') karta
1657 card=card(:ilen(card)+1)//karta
1661 subroutine read_dist_constr
1662 implicit real*8 (a-h,o-z)
1663 include 'DIMENSIONS'
1667 include 'COMMON.SETUP'
1668 include 'COMMON.CONTROL'
1669 include 'COMMON.CHAIN'
1670 include 'COMMON.IOUNITS'
1671 include 'COMMON.SBRIDGE'
1672 integer ifrag_(2,100),ipair_(2,100)
1673 double precision wfrag_(100),wpair_(100)
1674 character*500 controlcard
1675 write (iout,*) "Calling read_dist_constr"
1676 write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1678 call card_concat(controlcard)
1679 call readi(controlcard,"NFRAG",nfrag_,0)
1680 call readi(controlcard,"NPAIR",npair_,0)
1681 call readi(controlcard,"NDIST",ndist_,0)
1682 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1683 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1684 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1685 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1686 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1687 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1688 c write (iout,*) "IFRAG"
1690 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1692 c write (iout,*) "IPAIR"
1694 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1698 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1699 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1700 & ifrag_(2,i)=nstart_sup+nsup-1
1701 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1703 if (wfrag_(i).gt.0.0d0) then
1704 do j=ifrag_(1,i),ifrag_(2,i)-1
1705 do k=j+1,ifrag_(2,i)
1706 write (iout,*) "j",j," k",k
1708 if (constr_dist.eq.1) then
1713 forcon(nhpb)=wfrag_(i)
1714 else if (constr_dist.eq.2) then
1715 if (ddjk.le.dist_cut) then
1720 forcon(nhpb)=wfrag_(i)
1727 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1730 if (.not.out1file .or. me.eq.king)
1731 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1732 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1734 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1735 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1742 if (wpair_(i).gt.0.0d0) then
1750 do j=ifrag_(1,ii),ifrag_(2,ii)
1751 do k=ifrag_(1,jj),ifrag_(2,jj)
1755 forcon(nhpb)=wpair_(i)
1756 dhpb(nhpb)=dist(j,k)
1758 if (.not.out1file .or. me.eq.king)
1759 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1760 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1762 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1763 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1770 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1771 if (forcon(nhpb+1).gt.0.0d0) then
1773 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1775 if (.not.out1file .or. me.eq.king)
1776 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1777 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1779 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1780 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1787 c-------------------------------------------------------------------------------
1789 subroutine flush(iu)
1794 subroutine flush(iu)
1799 c------------------------------------------------------------------------------
1800 subroutine copy_to_tmp(source)
1801 include "DIMENSIONS"
1802 include "COMMON.IOUNITS"
1803 character*(*) source
1804 character* 256 tmpfile
1808 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1809 inquire(file=tmpfile,exist=ex)
1811 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1812 & " to temporary directory..."
1813 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1814 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1818 c------------------------------------------------------------------------------
1819 subroutine move_from_tmp(source)
1820 include "DIMENSIONS"
1821 include "COMMON.IOUNITS"
1822 character*(*) source
1825 write (*,*) "Moving ",source(:ilen(source)),
1826 & " from temporary directory to working directory"
1827 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1828 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1831 c------------------------------------------------------------------------------
1832 subroutine random_init(seed)
1834 C Initialize random number generator
1836 implicit real*8 (a-h,o-z)
1837 include 'DIMENSIONS'
1843 logical OKRandom, prng_restart
1845 integer iseed_array(4)
1847 include 'COMMON.IOUNITS'
1848 include 'COMMON.TIME1'
1849 c include 'COMMON.THREAD'
1850 include 'COMMON.SBRIDGE'
1851 include 'COMMON.CONTROL'
1852 include 'COMMON.MCM'
1853 c include 'COMMON.MAP'
1854 include 'COMMON.HEADER'
1855 c include 'COMMON.CSA'
1856 include 'COMMON.CHAIN'
1857 c include 'COMMON.MUCA'
1858 c include 'COMMON.MD'
1859 include 'COMMON.FFIELD'
1860 include 'COMMON.SETUP'
1861 iseed=-dint(dabs(seed))
1862 if (iseed.eq.0) then
1863 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1864 & 'Random seed undefined. The program will stop.'
1865 write (*,'(/80(1h*)/20x,a/80(1h*))')
1866 & 'Random seed undefined. The program will stop.'
1868 call mpi_finalize(mpi_comm_world,ierr)
1870 stop 'Bad random seed.'
1873 if (fg_rank.eq.0) then
1877 if(me.eq.king .or. .not. out1file)
1878 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1879 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1880 OKRandom = prng_restart(me,iseedi8)
1883 tmp=65536.0d0**(4-i)
1884 iseed_array(i) = dint(seed/tmp)
1885 seed=seed-iseed_array(i)*tmp
1887 if(me.eq.king .or. .not. out1file)
1888 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1889 & (iseed_array(i),i=1,4)
1890 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1891 & (iseed_array(i),i=1,4)
1892 OKRandom = prng_restart(me,iseed_array)
1895 r1=ran_number(0.0D0,1.0D0)
1896 if(me.eq.king .or. .not. out1file)
1897 & write (iout,*) 'ran_num',r1
1898 if (r1.lt.0.0d0) OKRandom=.false.
1900 if (.not.OKRandom) then
1901 write (iout,*) 'PRNG IS NOT WORKING!!!'
1902 print *,'PRNG IS NOT WORKING!!!'
1905 call mpi_abort(mpi_comm_world,error_msg,ierr)
1908 write (iout,*) 'too many processors for parallel prng'
1909 write (*,*) 'too many processors for parallel prng'
1917 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)