2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
16 C Read control parameters for energy minimzation if required
17 if (minim) call read_minim
18 C Read MCM control parameters if required
19 c if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 c if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 c if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 c if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 if (modecalc.eq.8) then
32 inquire (file="fort.40",exist=file_exist)
33 if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 c include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
74 c include 'COMMON.MAP'
75 include 'COMMON.HEADER'
76 c include 'COMMON.CSA'
77 include 'COMMON.CHAIN'
78 c include 'COMMON.MUCA'
80 include 'COMMON.FFIELD'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
100 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
101 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
102 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
103 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
104 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
105 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
106 call reada(controlcard,'DRMS',drms,0.1D0)
107 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
108 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
109 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
110 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
111 write (iout,'(a,f10.1)')'DRMS = ',drms
112 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
113 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
115 call readi(controlcard,'NZ_START',nz_start,0)
116 call readi(controlcard,'NZ_END',nz_end,0)
117 call readi(controlcard,'IZ_SC',iz_sc,0)
119 safety = 60.0d0*safety
122 call reada(controlcard,"T_BATH",t_bath,300.0d0)
123 minim=(index(controlcard,'MINIMIZE').gt.0)
124 dccart=(index(controlcard,'CART').gt.0)
125 overlapsc=(index(controlcard,'OVERLAP').gt.0)
126 overlapsc=.not.overlapsc
127 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
128 searchsc=.not.searchsc
129 sideadd=(index(controlcard,'SIDEADD').gt.0)
130 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
131 outpdb=(index(controlcard,'PDBOUT').gt.0)
132 outmol2=(index(controlcard,'MOL2OUT').gt.0)
133 pdbref=(index(controlcard,'PDBREF').gt.0)
134 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
135 indpdb=index(controlcard,'PDBSTART')
136 extconf=(index(controlcard,'EXTCONF').gt.0)
137 call readi(controlcard,'IPRINT',iprint,0)
138 call readi(controlcard,'MAXGEN',maxgen,10000)
139 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
140 call readi(controlcard,"KDIAG",kdiag,0)
141 call readi(controlcard,"RESCALE_MODE",rescale_mode,1)
142 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
143 & write (iout,*) "RESCALE_MODE",rescale_mode
144 split_ene=index(controlcard,'SPLIT_ENE').gt.0
145 if (index(controlcard,'REGULAR').gt.0.0D0) then
146 call reada(controlcard,'WEIDIS',weidis,0.1D0)
150 if (index(controlcard,'CHECKGRAD').gt.0) then
152 if (index(controlcard,'CART').gt.0) then
154 elseif (index(controlcard,'CARINT').gt.0) then
159 elseif (index(controlcard,'THREAD').gt.0) then
161 call readi(controlcard,'THREAD',nthread,0)
162 if (nthread.gt.0) then
163 call reada(controlcard,'WEIDIS',weidis,0.1D0)
166 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
167 stop 'Error termination in Read_Control.'
169 else if (index(controlcard,'MCMA').gt.0) then
171 else if (index(controlcard,'MCEE').gt.0) then
173 else if (index(controlcard,'MULTCONF').gt.0) then
175 else if (index(controlcard,'MAP').gt.0) then
177 call readi(controlcard,'MAP',nmap,0)
178 else if (index(controlcard,'CSA').gt.0) then
180 crc else if (index(controlcard,'ZSCORE').gt.0) then
182 crc ZSCORE is rm from UNRES, modecalc=9 is available
185 cfcm else if (index(controlcard,'MCMF').gt.0) then
187 else if (index(controlcard,'SOFTREG').gt.0) then
189 else if (index(controlcard,'CHECK_BOND').gt.0) then
191 else if (index(controlcard,'TEST').gt.0) then
193 else if (index(controlcard,'MD').gt.0) then
195 else if (index(controlcard,'RE ').gt.0) then
199 lmuca=index(controlcard,'MUCA').gt.0
200 call readi(controlcard,'MUCADYN',mucadyn,0)
201 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
202 if (lmuca .and. (me.eq.king .or. .not.out1file ))
204 write (iout,*) 'MUCADYN=',mucadyn
205 write (iout,*) 'MUCASMOOTH=',muca_smooth
208 iscode=index(controlcard,'ONE_LETTER')
209 indphi=index(controlcard,'PHI')
210 indback=index(controlcard,'BACK')
211 iranconf=index(controlcard,'RAND_CONF')
212 i2ndstr=index(controlcard,'USE_SEC_PRED')
213 gradout=index(controlcard,'GRADOUT').gt.0
214 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 if(me.eq.king.or..not.out1file)
217 & write (iout,'(2a)') diagmeth(kdiag),
218 & ' routine used to diagonalize matrices.'
221 c------------------------------------------------------------------------------
224 C Read molecular data.
226 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
235 include 'COMMON.INTERACT'
236 include 'COMMON.LOCAL'
237 include 'COMMON.NAMES'
238 include 'COMMON.CHAIN'
239 include 'COMMON.FFIELD'
240 include 'COMMON.SBRIDGE'
241 include 'COMMON.HEADER'
242 include 'COMMON.CONTROL'
243 c include 'COMMON.DBASE'
244 c include 'COMMON.THREAD'
245 include 'COMMON.CONTACTS'
246 include 'COMMON.TORCNSTR'
247 include 'COMMON.TIME1'
248 include 'COMMON.BOUNDS'
249 c include 'COMMON.MD'
250 c include 'COMMON.REMD'
251 include 'COMMON.SETUP'
252 character*4 sequence(maxres)
254 double precision x(maxvar)
255 character*256 pdbfile
256 character*320 weightcard
257 character*80 weightcard_t,ucase
258 dimension itype_pdb(maxres)
259 common /pizda/ itype_pdb
260 logical seq_comp,fail
261 double precision energia(0:n_ene)
267 C Read weights of the subsequent energy terms.
269 call card_concat(weightcard)
270 call reada(weightcard,'WLONG',wlong,1.0D0)
271 call reada(weightcard,'WSC',wsc,wlong)
272 call reada(weightcard,'WSCP',wscp,wlong)
273 call reada(weightcard,'WELEC',welec,1.0D0)
274 call reada(weightcard,'WVDWPP',wvdwpp,welec)
275 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
276 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
277 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
278 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
279 call reada(weightcard,'WTURN3',wturn3,1.0D0)
280 call reada(weightcard,'WTURN4',wturn4,1.0D0)
281 call reada(weightcard,'WTURN6',wturn6,1.0D0)
282 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
283 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
284 call reada(weightcard,'WBOND',wbond,1.0D0)
285 call reada(weightcard,'WTOR',wtor,1.0D0)
286 call reada(weightcard,'WTORD',wtor_d,1.0D0)
287 call reada(weightcard,'WANG',wang,1.0D0)
288 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
290 call reada(weightcard,'WDFAD',wdfa_dist,1.0d0)
291 call reada(weightcard,'WDFAT',wdfa_tor,1.0d0)
292 call reada(weightcard,'WDFAN',wdfa_nei,1.0d0)
293 call reada(weightcard,'WDFAB',wdfa_beta,1.0d0)
295 call reada(weightcard,'SCAL14',scal14,0.4D0)
296 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
297 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
298 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
299 call reada(weightcard,'TEMP0',temp0,300.0d0)
300 if (index(weightcard,'SOFT').gt.0) ipot=6
301 C 12/1/95 Added weight for the multi-body term WCORR
302 call reada(weightcard,'WCORRH',wcorr,1.0D0)
303 if (wcorr4.gt.0.0d0) wcorr=wcorr4
324 weights(24)=wdfa_dist
327 weights(27)=wdfa_beta
330 if(me.eq.king.or..not.out1file)
331 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
332 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
334 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
336 10 format (/'Energy-term weights (unscaled):'//
337 & 'WSCC= ',f10.6,' (SC-SC)'/
338 & 'WSCP= ',f10.6,' (SC-p)'/
339 & 'WELEC= ',f10.6,' (p-p electr)'/
340 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
341 & 'WBOND= ',f10.6,' (stretching)'/
342 & 'WANG= ',f10.6,' (bending)'/
343 & 'WSCLOC= ',f10.6,' (SC local)'/
344 & 'WTOR= ',f10.6,' (torsional)'/
345 & 'WTORD= ',f10.6,' (double torsional)'/
346 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
347 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
348 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
349 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
350 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
351 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
352 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
353 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
354 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
355 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
356 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
357 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
358 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
360 if(me.eq.king.or..not.out1file)then
361 if (wcorr4.gt.0.0d0) then
362 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
363 & 'between contact pairs of peptide groups'
364 write (iout,'(2(a,f5.3/))')
365 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
366 & 'Range of quenching the correlation terms:',2*delt_corr
367 else if (wcorr.gt.0.0d0) then
368 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
369 & 'between contact pairs of peptide groups'
371 write (iout,'(a,f8.3)')
372 & 'Scaling factor of 1,4 SC-p interactions:',scal14
373 write (iout,'(a,f8.3)')
374 & 'General scaling factor of SC-p interactions:',scalscp
376 r0_corr=cutoff_corr-delt_corr
378 aad(i,1)=scalscp*aad(i,1)
379 aad(i,2)=scalscp*aad(i,2)
380 bad(i,1)=scalscp*bad(i,1)
381 bad(i,2)=scalscp*bad(i,2)
383 call rescale_weights(t_bath)
384 if(me.eq.king.or..not.out1file)
385 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
386 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
388 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
390 22 format (/'Energy-term weights (scaled):'//
391 & 'WSCC= ',f10.6,' (SC-SC)'/
392 & 'WSCP= ',f10.6,' (SC-p)'/
393 & 'WELEC= ',f10.6,' (p-p electr)'/
394 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
395 & 'WBOND= ',f10.6,' (stretching)'/
396 & 'WANG= ',f10.6,' (bending)'/
397 & 'WSCLOC= ',f10.6,' (SC local)'/
398 & 'WTOR= ',f10.6,' (torsional)'/
399 & 'WTORD= ',f10.6,' (double torsional)'/
400 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
401 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
402 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
403 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
404 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
405 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
406 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
407 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
408 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
409 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
410 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
411 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
412 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
414 if(me.eq.king.or..not.out1file)
415 & write (iout,*) "Reference temperature for weights calculation:",
417 call reada(weightcard,"D0CM",d0cm,3.78d0)
418 call reada(weightcard,"AKCM",akcm,15.1d0)
419 call reada(weightcard,"AKTH",akth,11.0d0)
420 call reada(weightcard,"AKCT",akct,12.0d0)
421 call reada(weightcard,"V1SS",v1ss,-1.08d0)
422 call reada(weightcard,"V2SS",v2ss,7.61d0)
423 call reada(weightcard,"V3SS",v3ss,13.7d0)
424 call reada(weightcard,"EBR",ebr,-5.50D0)
425 if(me.eq.king.or..not.out1file) then
426 write (iout,*) "Parameters of the SS-bond potential:"
427 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
429 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
430 write (iout,*) "EBR",ebr
431 print *,'indpdb=',indpdb,' pdbref=',pdbref
433 if (indpdb.gt.0 .or. pdbref) then
434 read(inp,'(a)') pdbfile
435 if(me.eq.king.or..not.out1file)
436 & write (iout,'(2a)') 'PDB data will be read from file ',
437 & pdbfile(:ilen(pdbfile))
438 open(ipdbin,file=pdbfile,status='old',err=33)
440 33 write (iout,'(a)') 'Error opening PDB file.'
443 c print *,'Begin reading pdb data'
445 c print *,'Finished reading pdb data'
446 if(me.eq.king.or..not.out1file)
447 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
448 & ' nstart_sup=',nstart_sup
450 itype_pdb(i)=itype(i)
454 nct=nstart_sup+nsup-1
455 call contact(.false.,ncont_ref,icont_ref,co)
458 if(me.eq.king.or..not.out1file)
459 & write(iout,*)'Adding sidechains'
466 do while (fail.and.nsi.le.maxsi)
467 c call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
470 if(fail) write(iout,*)'Adding sidechain failed for res ',
471 & i,' after ',nsi,' trials'
477 if (indpdb.eq.0) then
478 C Read sequence if not taken from the pdb file.
480 c print *,'nres=',nres
481 if (iscode.gt.0) then
482 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
484 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
486 C Convert sequence to numeric code
488 itype(i)=rescode(i,sequence(i),iscode)
490 C Assign initial virtual bond lengths
496 vbld(i+nres)=dsc(itype(i))
497 vbld_inv(i+nres)=dsc_inv(itype(i))
498 c write (iout,*) "i",i," itype",itype(i),
499 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
503 c print '(20i4)',(itype(i),i=1,nres)
506 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
508 if (itype(i).eq.21) then
512 else if (itype(i+1).ne.20) then
514 else if (itype(i).ne.20) then
521 if(me.eq.king.or..not.out1file)then
522 write (iout,*) "ITEL"
524 write (iout,*) i,itype(i),itel(i)
526 print *,'Call Read_Bridge.'
529 C 8/13/98 Set limits to generating the dihedral angles
534 read (inp,*) ndih_constr
535 if (ndih_constr.gt.0) then
537 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
538 if(me.eq.king.or..not.out1file)then
540 & 'There are',ndih_constr,' constraints on phi angles.'
542 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
546 phi0(i)=deg2rad*phi0(i)
547 drange(i)=deg2rad*drange(i)
549 if(me.eq.king.or..not.out1file)
550 & write (iout,*) 'FTORS',ftors
553 phibound(1,ii) = phi0(i)-drange(i)
554 phibound(2,ii) = phi0(i)+drange(i)
561 write (iout,'(a)') 'Boundaries in phi angle sampling:'
563 write (iout,'(a3,i5,2f10.1)')
564 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
570 cd print *,'NNT=',NNT,' NCT=',NCT
571 if (itype(1).eq.21) nnt=2
572 if (itype(nres).eq.21) nct=nct-1
574 C Juyong:READ init_vars
575 C Initialize variables!
576 C Juyong:READ read_info
577 C READ fragment information!!
578 C both routines should be in dfa.F file!!
581 print*, 'init_dfa_vars finished!'
583 print*, 'read_dfa_info finished!'
589 if(me.eq.king.or..not.out1file)
590 & write (iout,'(a,i3)') 'nsup=',nsup
592 if (nsup.le.(nct-nnt+1)) then
593 do i=0,nct-nnt+1-nsup
594 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
600 & 'Error - sequences to be superposed do not match.'
603 do i=0,nsup-(nct-nnt+1)
604 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
606 nstart_sup=nstart_sup+i
612 & 'Error - sequences to be superposed do not match.'
615 if (nsup.eq.0) nsup=nct-nnt
616 if (nstart_sup.eq.0) nstart_sup=nnt
617 if (nstart_seq.eq.0) nstart_seq=nnt
618 if(me.eq.king.or..not.out1file)
619 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
620 & ' nstart_seq=',nstart_seq
622 c--- Zscore rms -------
623 if (nz_start.eq.0) nz_start=nnt
624 if (nz_end.eq.0 .and. nsup.gt.0) then
626 else if (nz_end.eq.0) then
629 if(me.eq.king.or..not.out1file)then
630 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
631 write (iout,*) 'IZ_SC=',iz_sc
633 c----------------------
636 if (.not.pdbref) then
637 call read_angles(inp,*38)
639 38 write (iout,'(a)') 'Error reading reference structure.'
641 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
642 stop 'Error reading reference structure'
646 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
655 call contact(.true.,ncont_ref,icont_ref,co)
657 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
659 if (constr_dist.gt.0) call read_dist_constr
660 c write (iout,*) "After read_dist_constr nhpb",nhpb
662 if(me.eq.king.or..not.out1file)
663 & write (iout,*) 'Contact order:',co
665 if(me.eq.king.or..not.out1file)
666 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
669 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
671 if(me.eq.king.or..not.out1file)
672 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
673 & icont_ref(1,i),' ',
674 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
678 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
679 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
680 & modecalc.ne.10) then
681 C If input structure hasn't been supplied from the PDB file read or generate
683 if (iranconf.eq.0 .and. .not. extconf) then
684 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
685 & write (iout,'(a)') 'Initial geometry will be read in.'
687 read(inp,'(8f10.5)',end=36,err=36)
688 & ((c(l,k),l=1,3),k=1,nres),
689 & ((c(l,k+nres),l=1,3),k=nnt,nct)
690 call int_from_cart1(.false.)
693 dc(j,i)=c(j,i+1)-c(j,i)
694 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
698 if (itype(i).ne.10) then
700 dc(j,i+nres)=c(j,i+nres)-c(j,i)
701 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
707 call read_angles(inp,*36)
710 36 write (iout,'(a)') 'Error reading angle file.'
712 call mpi_finalize( MPI_COMM_WORLD,IERR )
714 stop 'Error reading angle file.'
716 else if (extconf) then
717 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
718 & write (iout,'(a)') 'Extended chain initial geometry.'
720 theta(i)=90d0*deg2rad
726 alph(i)=110d0*deg2rad
729 omeg(i)=-120d0*deg2rad
732 if(me.eq.king.or..not.out1file)
733 & write (iout,'(a)') 'Random-generated initial geometry.'
737 if (me.eq.king .or. fg_rank.eq.0 .and. (
738 & modecalc.eq.12 .or. modecalc.eq.14) ) then
742 c call gen_rand_conf(itmp,*30)
744 30 write (iout,*) 'Failed to generate random conformation',
746 write (*,*) 'Processor:',me,
747 & ' Failed to generate random conformation',
757 write (iout,'(a,i3,a)') 'Processor:',me,
758 & ' error in generating random conformation.'
759 write (*,'(a,i3,a)') 'Processor:',me,
760 & ' error in generating random conformation.'
763 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
769 c call gen_rand_conf(itmp,*31)
771 31 write (iout,*) 'Failed to generate random conformation',
773 write (*,*) 'Failed to generate random conformation',
776 write (iout,'(a,i3,a)') 'Processor:',me,
777 & ' error in generating random conformation.'
778 write (*,'(a,i3,a)') 'Processor:',me,
779 & ' error in generating random conformation.'
784 elseif (modecalc.eq.4) then
785 read (inp,'(a)') intinname
786 open (intin,file=intinname,status='old',err=333)
787 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
788 & write (iout,'(a)') 'intinname',intinname
789 write (*,'(a)') 'Processor',myrank,' intinname',intinname
791 333 write (iout,'(2a)') 'Error opening angle file ',intinname
793 call MPI_Finalize(MPI_COMM_WORLD,IERR)
795 stop 'Error opening angle file.'
799 C Generate distance constraints, if the PDB structure is to be regularized.
800 if (nthread.gt.0) then
804 if (me.eq.king .or. .not. out1file)
806 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
807 write (iout,'(/a,i3,a)')
808 & 'The chain contains',ns,' disulfide-bridging cysteines.'
809 write (iout,'(20i4)') (iss(i),i=1,ns)
810 write (iout,'(/a/)') 'Pre-formed links are:'
816 if (me.eq.king.or..not.out1file)
817 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
818 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
823 c if (i2ndstr.gt.0) call secstrp2dihc
824 c call geom_to_var(nvar,x)
825 c call etotal(energia(0))
826 c call enerprint(energia(0))
827 c call briefout(0,etot)
829 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
830 cd write (iout,'(a)') 'Variable list:'
831 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
833 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
834 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
835 & 'Processor',myrank,': end reading molecular data.'
839 c--------------------------------------------------------------------------
840 logical function seq_comp(itypea,itypeb,length)
842 integer length,itypea(length),itypeb(length)
845 if (itypea(i).ne.itypeb(i)) then
853 c-----------------------------------------------------------------------------
854 subroutine read_bridge
855 C Read information about disulfide bridges.
856 implicit real*8 (a-h,o-z)
861 include 'COMMON.IOUNITS'
864 include 'COMMON.INTERACT'
865 include 'COMMON.LOCAL'
866 include 'COMMON.NAMES'
867 include 'COMMON.CHAIN'
868 include 'COMMON.FFIELD'
869 include 'COMMON.SBRIDGE'
870 include 'COMMON.HEADER'
871 include 'COMMON.CONTROL'
872 c include 'COMMON.DBASE'
873 c include 'COMMON.THREAD'
874 include 'COMMON.TIME1'
875 include 'COMMON.SETUP'
876 C Read bridging residues.
877 read (inp,*) ns,(iss(i),i=1,ns)
879 if(me.eq.king.or..not.out1file)
880 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
881 C Check whether the specified bridging residues are cystines.
883 if (itype(iss(i)).ne.1) then
884 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
885 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
886 & ' can form a disulfide bridge?!!!'
887 write (*,'(2a,i3,a)')
888 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
889 & ' can form a disulfide bridge?!!!'
891 call MPI_Finalize(MPI_COMM_WORLD,ierror)
896 C Read preformed bridges.
898 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
899 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
902 C Check if the residues involved in bridges are in the specified list of
906 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
907 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
908 write (iout,'(a,i3,a)') 'Disulfide pair',i,
909 & ' contains residues present in other pairs.'
910 write (*,'(a,i3,a)') 'Disulfide pair',i,
911 & ' contains residues present in other pairs.'
913 call MPI_Finalize(MPI_COMM_WORLD,ierror)
919 if (ihpb(i).eq.iss(j)) goto 10
921 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
924 if (jhpb(i).eq.iss(j)) goto 20
926 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
939 c----------------------------------------------------------------------------
940 subroutine read_x(kanal,*)
941 implicit real*8 (a-h,o-z)
945 include 'COMMON.CHAIN'
946 include 'COMMON.IOUNITS'
947 include 'COMMON.CONTROL'
948 include 'COMMON.LOCAL'
949 include 'COMMON.INTERACT'
950 c Read coordinates from input
952 read(kanal,'(8f10.5)',end=10,err=10)
953 & ((c(l,k),l=1,3),k=1,nres),
954 & ((c(l,k+nres),l=1,3),k=nnt,nct)
957 c(j,2*nres)=c(j,nres)
959 call int_from_cart1(.false.)
962 dc(j,i)=c(j,i+1)-c(j,i)
963 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
967 if (itype(i).ne.10) then
969 dc(j,i+nres)=c(j,i+nres)-c(j,i)
970 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
978 c------------------------------------------------------------------------------
980 implicit real*8 (a-h,o-z)
982 include 'COMMON.IOUNITS'
985 include 'COMMON.INTERACT'
986 include 'COMMON.LOCAL'
987 include 'COMMON.NAMES'
988 include 'COMMON.CHAIN'
989 include 'COMMON.FFIELD'
990 include 'COMMON.SBRIDGE'
991 include 'COMMON.HEADER'
992 include 'COMMON.CONTROL'
993 c include 'COMMON.DBASE'
994 c include 'COMMON.THREAD'
995 include 'COMMON.TIME1'
996 C Set up variable list.
1002 if (itype(i).ne.10) then
1004 ialph(i,1)=nvar+nside
1008 if (indphi.gt.0) then
1010 else if (indback.gt.0) then
1015 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1018 c----------------------------------------------------------------------------
1019 subroutine gen_dist_constr
1020 C Generate CA distance constraints.
1021 implicit real*8 (a-h,o-z)
1022 include 'DIMENSIONS'
1023 include 'COMMON.IOUNITS'
1024 include 'COMMON.GEO'
1025 include 'COMMON.VAR'
1026 include 'COMMON.INTERACT'
1027 include 'COMMON.LOCAL'
1028 include 'COMMON.NAMES'
1029 include 'COMMON.CHAIN'
1030 include 'COMMON.FFIELD'
1031 include 'COMMON.SBRIDGE'
1032 include 'COMMON.HEADER'
1033 include 'COMMON.CONTROL'
1034 c include 'COMMON.DBASE'
1035 c include 'COMMON.THREAD'
1036 include 'COMMON.TIME1'
1037 dimension itype_pdb(maxres)
1038 common /pizda/ itype_pdb
1040 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1041 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1042 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1044 do i=nstart_sup,nstart_sup+nsup-1
1045 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1046 cd & ' seq_pdb', restyp(itype_pdb(i))
1047 do j=i+2,nstart_sup+nsup-1
1049 ihpb(nhpb)=i+nstart_seq-nstart_sup
1050 jhpb(nhpb)=j+nstart_seq-nstart_sup
1052 dhpb(nhpb)=dist(i,j)
1055 cd write (iout,'(a)') 'Distance constraints:'
1060 cd if (ii.gt.nres) then
1065 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1066 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1067 cd & dhpb(i),forcon(i)
1071 c----------------------------------------------------------------------------
1073 implicit real*8 (a-h,o-z)
1074 include 'DIMENSIONS'
1075 include 'COMMON.IOUNITS'
1076 include 'COMMON.GEO'
1077 include 'COMMON.CSA'
1078 include 'COMMON.BANK'
1079 include 'COMMON.CONTROL'
1081 character*620 mcmcard
1082 call card_concat(mcmcard)
1084 call readi(mcmcard,'NCONF',nconf,50)
1085 call readi(mcmcard,'NADD',nadd,0)
1086 call readi(mcmcard,'JSTART',jstart,1)
1087 call readi(mcmcard,'JEND',jend,1)
1088 call readi(mcmcard,'NSTMAX',nstmax,500000)
1089 call readi(mcmcard,'N0',n0,1)
1090 call readi(mcmcard,'N1',n1,6)
1091 call readi(mcmcard,'N2',n2,4)
1092 call readi(mcmcard,'N3',n3,0)
1093 call readi(mcmcard,'N4',n4,0)
1094 call readi(mcmcard,'N5',n5,0)
1095 call readi(mcmcard,'N6',n6,10)
1096 call readi(mcmcard,'N7',n7,0)
1097 call readi(mcmcard,'N8',n8,0)
1098 call readi(mcmcard,'N9',n9,0)
1099 call readi(mcmcard,'N14',n14,0)
1100 call readi(mcmcard,'N15',n15,0)
1101 call readi(mcmcard,'N16',n16,0)
1102 call readi(mcmcard,'N17',n17,0)
1103 call readi(mcmcard,'N18',n18,0)
1105 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1107 call readi(mcmcard,'NDIFF',ndiff,2)
1108 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1109 call readi(mcmcard,'IS1',is1,1)
1110 call readi(mcmcard,'IS2',is2,8)
1111 call readi(mcmcard,'NRAN0',nran0,4)
1112 call readi(mcmcard,'NRAN1',nran1,2)
1113 call readi(mcmcard,'IRR',irr,1)
1114 call readi(mcmcard,'NSEED',nseed,20)
1115 call readi(mcmcard,'NTOTAL',ntotal,10000)
1116 call reada(mcmcard,'CUT1',cut1,2.0d0)
1117 call reada(mcmcard,'CUT2',cut2,5.0d0)
1118 call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1119 call readi(mcmcard,'ICMAX',icmax,1)
1120 call readi(mcmcard,'NBANKM',nbankm,400)
1121 call readi(mcmcard,'IUCUT',iucut,2)
1122 call readi(mcmcard,'IRESTART',irestart,0)
1123 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1126 call reada(mcmcard,'DELE',dele,20.0d0)
1127 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1128 call readi(mcmcard,'IREF',iref,0)
1129 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1130 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1131 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1132 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1133 write (iout,*) "NCONF_IN",nconf_in
1137 subroutine read_minim
1138 implicit real*8 (a-h,o-z)
1139 include 'DIMENSIONS'
1140 include 'COMMON.MINIM'
1141 include 'COMMON.IOUNITS'
1143 character*320 minimcard
1144 call card_concat(minimcard)
1145 call readi(minimcard,'MAXMIN',maxmin,2000)
1146 call readi(minimcard,'MAXFUN',maxfun,5000)
1147 call readi(minimcard,'MINMIN',minmin,maxmin)
1148 call readi(minimcard,'MINFUN',minfun,maxmin)
1149 call reada(minimcard,'TOLF',tolf,1.0D-2)
1150 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1151 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1152 & 'Options in energy minimization:'
1153 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1154 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1155 & 'MinMin:',MinMin,' MinFun:',MinFun,
1156 & ' TolF:',TolF,' RTolF:',RTolF
1159 c----------------------------------------------------------------------------
1160 subroutine read_angles(kanal,*)
1161 implicit real*8 (a-h,o-z)
1162 include 'DIMENSIONS'
1163 include 'COMMON.GEO'
1164 include 'COMMON.VAR'
1165 include 'COMMON.CHAIN'
1166 include 'COMMON.IOUNITS'
1167 include 'COMMON.CONTROL'
1168 c Read angles from input
1170 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1171 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1172 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1173 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1176 c 9/7/01 avoid 180 deg valence angle
1177 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1179 theta(i)=deg2rad*theta(i)
1180 phi(i)=deg2rad*phi(i)
1181 alph(i)=deg2rad*alph(i)
1182 omeg(i)=deg2rad*omeg(i)
1187 c----------------------------------------------------------------------------
1188 subroutine reada(rekord,lancuch,wartosc,default)
1190 character*(*) rekord,lancuch
1191 double precision wartosc,default
1194 iread=index(rekord,lancuch)
1195 if (iread.eq.0) then
1199 iread=iread+ilen(lancuch)+1
1200 read (rekord(iread:),*,err=10,end=10) wartosc
1205 c----------------------------------------------------------------------------
1206 subroutine readi(rekord,lancuch,wartosc,default)
1208 character*(*) rekord,lancuch
1209 integer 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 multreadi(rekord,lancuch,tablica,dim,default)
1227 integer tablica(dim),default
1228 character*(*) rekord,lancuch
1235 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1236 if (iread.eq.0) return
1237 iread=iread+ilen(lancuch)+1
1238 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1241 c----------------------------------------------------------------------------
1242 subroutine multreada(rekord,lancuch,tablica,dim,default)
1245 double precision 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 openunits
1261 implicit real*8 (a-h,o-z)
1262 include 'DIMENSIONS'
1265 character*16 form,nodename
1268 include 'COMMON.SETUP'
1269 include 'COMMON.IOUNITS'
1270 c include 'COMMON.MD'
1271 include 'COMMON.CONTROL'
1272 integer lenpre,lenpot,ilen,lentmp
1274 character*3 out1file_text,ucase
1277 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1278 call getenv_loc("PREFIX",prefix)
1280 call getenv_loc("POT",pot)
1281 call getenv_loc("DIRTMP",tmpdir)
1282 call getenv_loc("CURDIR",curdir)
1283 call getenv_loc("OUT1FILE",out1file_text)
1284 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1285 out1file_text=ucase(out1file_text)
1286 if (out1file_text(1:1).eq."Y") then
1289 out1file=fg_rank.gt.0
1294 if (lentmp.gt.0) then
1295 write (*,'(80(1h!))')
1296 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1297 write (*,'(80(1h!))')
1298 write (*,*)"All output files will be on node /tmp directory."
1300 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1301 if (me.eq.king) then
1302 write (*,*) "The master node is ",nodename
1303 else if (fg_rank.eq.0) then
1304 write (*,*) "I am the CG slave node ",nodename
1306 write (*,*) "I am the FG slave node ",nodename
1309 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1310 lenpre = lentmp+lenpre+1
1312 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1313 C Get the names and open the input files
1314 #if defined(WINIFL) || defined(WINPGI)
1315 open(1,file=pref_orig(:ilen(pref_orig))//
1316 & '.inp',status='old',readonly,shared)
1317 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1318 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1319 C Get parameter filenames and open the parameter files.
1320 call getenv_loc('BONDPAR',bondname)
1321 open (ibond,file=bondname,status='old',readonly,shared)
1322 call getenv_loc('THETPAR',thetname)
1323 open (ithep,file=thetname,status='old',readonly,shared)
1325 call getenv_loc('THETPARPDB',thetname_pdb)
1326 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1328 call getenv_loc('ROTPAR',rotname)
1329 open (irotam,file=rotname,status='old',readonly,shared)
1331 call getenv_loc('ROTPARPDB',rotname_pdb)
1332 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1334 call getenv_loc('TORPAR',torname)
1335 open (itorp,file=torname,status='old',readonly,shared)
1336 call getenv_loc('TORDPAR',tordname)
1337 open (itordp,file=tordname,status='old',readonly,shared)
1338 call getenv_loc('FOURIER',fouriername)
1339 open (ifourier,file=fouriername,status='old',readonly,shared)
1340 call getenv_loc('ELEPAR',elename)
1341 open (ielep,file=elename,status='old',readonly,shared)
1342 call getenv_loc('SIDEPAR',sidename)
1343 open (isidep,file=sidename,status='old',readonly,shared)
1344 #elif (defined CRAY) || (defined AIX)
1345 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1347 c print *,"Processor",myrank," opened file 1"
1348 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1349 c print *,"Processor",myrank," opened file 9"
1350 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1351 C Get parameter filenames and open the parameter files.
1352 call getenv_loc('BONDPAR',bondname)
1353 open (ibond,file=bondname,status='old',action='read')
1354 c print *,"Processor",myrank," opened file IBOND"
1355 call getenv_loc('THETPAR',thetname)
1356 open (ithep,file=thetname,status='old',action='read')
1357 c print *,"Processor",myrank," opened file ITHEP"
1359 call getenv_loc('THETPARPDB',thetname_pdb)
1360 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1362 call getenv_loc('ROTPAR',rotname)
1363 open (irotam,file=rotname,status='old',action='read')
1364 c print *,"Processor",myrank," opened file IROTAM"
1366 call getenv_loc('ROTPARPDB',rotname_pdb)
1367 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1369 call getenv_loc('TORPAR',torname)
1370 open (itorp,file=torname,status='old',action='read')
1371 c print *,"Processor",myrank," opened file ITORP"
1372 call getenv_loc('TORDPAR',tordname)
1373 open (itordp,file=tordname,status='old',action='read')
1374 c print *,"Processor",myrank," opened file ITORDP"
1375 call getenv_loc('SCCORPAR',sccorname)
1376 open (isccor,file=sccorname,status='old',action='read')
1377 c print *,"Processor",myrank," opened file ISCCOR"
1378 call getenv_loc('FOURIER',fouriername)
1379 open (ifourier,file=fouriername,status='old',action='read')
1380 c print *,"Processor",myrank," opened file IFOURIER"
1381 call getenv_loc('ELEPAR',elename)
1382 open (ielep,file=elename,status='old',action='read')
1383 c print *,"Processor",myrank," opened file IELEP"
1384 call getenv_loc('SIDEPAR',sidename)
1385 open (isidep,file=sidename,status='old',action='read')
1386 c print *,"Processor",myrank," opened file ISIDEP"
1387 c print *,"Processor",myrank," opened parameter files"
1389 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1390 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1391 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1392 C Get parameter filenames and open the parameter files.
1393 call getenv_loc('BONDPAR',bondname)
1394 open (ibond,file=bondname,status='old')
1395 call getenv_loc('THETPAR',thetname)
1396 open (ithep,file=thetname,status='old')
1398 call getenv_loc('THETPARPDB',thetname_pdb)
1399 open (ithep_pdb,file=thetname_pdb,status='old')
1401 call getenv_loc('ROTPAR',rotname)
1402 open (irotam,file=rotname,status='old')
1404 call getenv_loc('ROTPARPDB',rotname_pdb)
1405 open (irotam_pdb,file=rotname_pdb,status='old')
1407 call getenv_loc('TORPAR',torname)
1408 open (itorp,file=torname,status='old')
1409 call getenv_loc('TORDPAR',tordname)
1410 open (itordp,file=tordname,status='old')
1411 call getenv_loc('SCCORPAR',sccorname)
1412 open (isccor,file=sccorname,status='old')
1413 call getenv_loc('FOURIER',fouriername)
1414 open (ifourier,file=fouriername,status='old')
1415 call getenv_loc('ELEPAR',elename)
1416 open (ielep,file=elename,status='old')
1417 call getenv_loc('SIDEPAR',sidename)
1418 open (isidep,file=sidename,status='old')
1420 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1422 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1423 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1424 C Get parameter filenames and open the parameter files.
1425 call getenv_loc('BONDPAR',bondname)
1426 open (ibond,file=bondname,status='old',readonly)
1427 call getenv_loc('THETPAR',thetname)
1428 open (ithep,file=thetname,status='old',readonly)
1430 call getenv_loc('THETPARPDB',thetname_pdb)
1431 print *,"thetname_pdb ",thetname_pdb
1432 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1433 print *,ithep_pdb," opened"
1435 call getenv_loc('ROTPAR',rotname)
1436 open (irotam,file=rotname,status='old',readonly)
1438 call getenv_loc('ROTPARPDB',rotname_pdb)
1439 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1441 call getenv_loc('TORPAR',torname)
1442 open (itorp,file=torname,status='old',readonly)
1443 call getenv_loc('TORDPAR',tordname)
1444 open (itordp,file=tordname,status='old',readonly)
1445 call getenv_loc('SCCORPAR',sccorname)
1446 open (isccor,file=sccorname,status='old',readonly)
1447 call getenv_loc('FOURIER',fouriername)
1448 open (ifourier,file=fouriername,status='old',readonly)
1449 call getenv_loc('ELEPAR',elename)
1450 open (ielep,file=elename,status='old',readonly)
1451 call getenv_loc('SIDEPAR',sidename)
1452 open (isidep,file=sidename,status='old',readonly)
1456 C 8/9/01 In the newest version SCp interaction constants are read from a file
1457 C Use -DOLDSCP to use hard-coded constants instead.
1459 call getenv_loc('SCPPAR',scpname)
1460 #if defined(WINIFL) || defined(WINPGI)
1461 open (iscpp,file=scpname,status='old',readonly,shared)
1462 #elif (defined CRAY) || (defined AIX)
1463 open (iscpp,file=scpname,status='old',action='read')
1465 open (iscpp,file=scpname,status='old')
1467 open (iscpp,file=scpname,status='old',readonly)
1470 call getenv_loc('PATTERN',patname)
1471 #if defined(WINIFL) || defined(WINPGI)
1472 open (icbase,file=patname,status='old',readonly,shared)
1473 #elif (defined CRAY) || (defined AIX)
1474 open (icbase,file=patname,status='old',action='read')
1476 open (icbase,file=patname,status='old')
1478 open (icbase,file=patname,status='old',readonly)
1481 C Open output file only for CG processes
1482 c print *,"Processor",myrank," fg_rank",fg_rank
1483 if (fg_rank.eq.0) then
1485 if (nodes.eq.1) then
1488 npos = dlog10(dfloat(nodes-1))+1
1490 if (npos.lt.3) npos=3
1491 write (liczba,'(i1)') npos
1492 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1494 write (liczba,form) me
1495 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1496 & liczba(:ilen(liczba))
1497 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1499 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1501 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1502 & liczba(:ilen(liczba))//'.mol2'
1503 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1504 & liczba(:ilen(liczba))//'.stat'
1506 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1507 & //liczba(:ilen(liczba))//'.stat')
1508 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1511 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1512 c & liczba(:ilen(liczba))//'.const'
1517 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1518 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1519 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1520 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1521 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1523 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1525 rest2name=prefix(:ilen(prefix))//'.rst'
1527 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1530 #if defined(AIX) || defined(PGI)
1531 if (me.eq.king .or. .not. out1file)
1532 & open(iout,file=outname,status='unknown')
1534 if (fg_rank.gt.0) then
1535 write (liczba,'(i3.3)') myrank/nfgtasks
1536 write (ll,'(bz,i3.3)') fg_rank
1537 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1542 open(igeom,file=intname,status='unknown',position='append')
1543 open(ipdb,file=pdbname,status='unknown')
1544 open(imol2,file=mol2name,status='unknown')
1545 open(istat,file=statname,status='unknown',position='append')
1547 c1out open(iout,file=outname,status='unknown')
1550 if (me.eq.king .or. .not.out1file)
1551 & open(iout,file=outname,status='unknown')
1553 if (fg_rank.gt.0) then
1554 write (liczba,'(i3.3)') myrank/nfgtasks
1555 write (ll,'(bz,i3.3)') fg_rank
1556 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1561 open(igeom,file=intname,status='unknown',access='append')
1562 open(ipdb,file=pdbname,status='unknown')
1563 open(imol2,file=mol2name,status='unknown')
1564 open(istat,file=statname,status='unknown',access='append')
1566 c1out open(iout,file=outname,status='unknown')
1569 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1570 csa_seed=prefix(:lenpre)//'.CSA.seed'
1571 csa_history=prefix(:lenpre)//'.CSA.history'
1572 csa_bank=prefix(:lenpre)//'.CSA.bank'
1573 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1574 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1575 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1576 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1577 csa_int=prefix(:lenpre)//'.int'
1578 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1579 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1580 csa_in=prefix(:lenpre)//'.CSA.in'
1581 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1584 write (iout,'(80(1h-))')
1585 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1586 write (iout,'(80(1h-))')
1587 write (iout,*) "Input file : ",
1588 & pref_orig(:ilen(pref_orig))//'.inp'
1589 write (iout,*) "Output file : ",
1590 & outname(:ilen(outname))
1592 write (iout,*) "Sidechain potential file : ",
1593 & sidename(:ilen(sidename))
1595 write (iout,*) "SCp potential file : ",
1596 & scpname(:ilen(scpname))
1598 write (iout,*) "Electrostatic potential file : ",
1599 & elename(:ilen(elename))
1600 write (iout,*) "Cumulant coefficient file : ",
1601 & fouriername(:ilen(fouriername))
1602 write (iout,*) "Torsional parameter file : ",
1603 & torname(:ilen(torname))
1604 write (iout,*) "Double torsional parameter file : ",
1605 & tordname(:ilen(tordname))
1606 write (iout,*) "SCCOR parameter file : ",
1607 & sccorname(:ilen(sccorname))
1608 write (iout,*) "Bond & inertia constant file : ",
1609 & bondname(:ilen(bondname))
1610 write (iout,*) "Bending parameter file : ",
1611 & thetname(:ilen(thetname))
1612 write (iout,*) "Rotamer parameter file : ",
1613 & rotname(:ilen(rotname))
1614 write (iout,*) "Threading database : ",
1615 & patname(:ilen(patname))
1617 &write (iout,*)" DIRTMP : ",
1619 write (iout,'(80(1h-))')
1623 c----------------------------------------------------------------------------
1624 subroutine card_concat(card)
1625 implicit real*8 (a-h,o-z)
1626 include 'DIMENSIONS'
1627 include 'COMMON.IOUNITS'
1629 character*80 karta,ucase
1631 read (inp,'(a)') karta
1634 do while (karta(80:80).eq.'&')
1635 card=card(:ilen(card)+1)//karta(:79)
1636 read (inp,'(a)') karta
1639 card=card(:ilen(card)+1)//karta
1643 subroutine read_dist_constr
1644 implicit real*8 (a-h,o-z)
1645 include 'DIMENSIONS'
1649 include 'COMMON.SETUP'
1650 include 'COMMON.CONTROL'
1651 include 'COMMON.CHAIN'
1652 include 'COMMON.IOUNITS'
1653 include 'COMMON.SBRIDGE'
1654 integer ifrag_(2,100),ipair_(2,100)
1655 double precision wfrag_(100),wpair_(100)
1656 character*500 controlcard
1657 write (iout,*) "Calling read_dist_constr"
1658 write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1660 call card_concat(controlcard)
1661 call readi(controlcard,"NFRAG",nfrag_,0)
1662 call readi(controlcard,"NPAIR",npair_,0)
1663 call readi(controlcard,"NDIST",ndist_,0)
1664 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1665 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1666 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1667 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1668 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1669 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1670 c write (iout,*) "IFRAG"
1672 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1674 c write (iout,*) "IPAIR"
1676 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1680 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1681 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1682 & ifrag_(2,i)=nstart_sup+nsup-1
1683 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1685 if (wfrag_(i).gt.0.0d0) then
1686 do j=ifrag_(1,i),ifrag_(2,i)-1
1687 do k=j+1,ifrag_(2,i)
1688 write (iout,*) "j",j," k",k
1690 if (constr_dist.eq.1) then
1695 forcon(nhpb)=wfrag_(i)
1696 else if (constr_dist.eq.2) then
1697 if (ddjk.le.dist_cut) then
1702 forcon(nhpb)=wfrag_(i)
1709 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1712 if (.not.out1file .or. me.eq.king)
1713 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1714 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1716 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1717 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1724 if (wpair_(i).gt.0.0d0) then
1732 do j=ifrag_(1,ii),ifrag_(2,ii)
1733 do k=ifrag_(1,jj),ifrag_(2,jj)
1737 forcon(nhpb)=wpair_(i)
1738 dhpb(nhpb)=dist(j,k)
1740 if (.not.out1file .or. me.eq.king)
1741 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1742 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1744 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1745 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1752 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1753 if (forcon(nhpb+1).gt.0.0d0) then
1755 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1757 if (.not.out1file .or. me.eq.king)
1758 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1759 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1761 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1762 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1769 c-------------------------------------------------------------------------------
1771 subroutine flush(iu)
1776 subroutine flush(iu)
1781 c------------------------------------------------------------------------------
1782 subroutine copy_to_tmp(source)
1783 include "DIMENSIONS"
1784 include "COMMON.IOUNITS"
1785 character*(*) source
1786 character* 256 tmpfile
1790 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1791 inquire(file=tmpfile,exist=ex)
1793 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1794 & " to temporary directory..."
1795 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1796 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1800 c------------------------------------------------------------------------------
1801 subroutine move_from_tmp(source)
1802 include "DIMENSIONS"
1803 include "COMMON.IOUNITS"
1804 character*(*) source
1807 write (*,*) "Moving ",source(:ilen(source)),
1808 & " from temporary directory to working directory"
1809 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1810 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1813 c------------------------------------------------------------------------------
1814 subroutine random_init(seed)
1816 C Initialize random number generator
1818 implicit real*8 (a-h,o-z)
1819 include 'DIMENSIONS'
1825 logical OKRandom, prng_restart
1827 integer iseed_array(4)
1829 include 'COMMON.IOUNITS'
1830 include 'COMMON.TIME1'
1831 c include 'COMMON.THREAD'
1832 include 'COMMON.SBRIDGE'
1833 include 'COMMON.CONTROL'
1834 include 'COMMON.MCM'
1835 c include 'COMMON.MAP'
1836 include 'COMMON.HEADER'
1837 c include 'COMMON.CSA'
1838 include 'COMMON.CHAIN'
1839 c include 'COMMON.MUCA'
1840 c include 'COMMON.MD'
1841 include 'COMMON.FFIELD'
1842 include 'COMMON.SETUP'
1843 iseed=-dint(dabs(seed))
1844 if (iseed.eq.0) then
1845 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1846 & 'Random seed undefined. The program will stop.'
1847 write (*,'(/80(1h*)/20x,a/80(1h*))')
1848 & 'Random seed undefined. The program will stop.'
1850 call mpi_finalize(mpi_comm_world,ierr)
1852 stop 'Bad random seed.'
1855 if (fg_rank.eq.0) then
1859 if(me.eq.king .or. .not. out1file)
1860 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1861 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1862 OKRandom = prng_restart(me,iseedi8)
1865 tmp=65536.0d0**(4-i)
1866 iseed_array(i) = dint(seed/tmp)
1867 seed=seed-iseed_array(i)*tmp
1869 if(me.eq.king .or. .not. out1file)
1870 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1871 & (iseed_array(i),i=1,4)
1872 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1873 & (iseed_array(i),i=1,4)
1874 OKRandom = prng_restart(me,iseed_array)
1877 r1=ran_number(0.0D0,1.0D0)
1878 if(me.eq.king .or. .not. out1file)
1879 & write (iout,*) 'ran_num',r1
1880 if (r1.lt.0.0d0) OKRandom=.false.
1882 if (.not.OKRandom) then
1883 write (iout,*) 'PRNG IS NOT WORKING!!!'
1884 print *,'PRNG IS NOT WORKING!!!'
1887 call mpi_abort(mpi_comm_world,error_msg,ierr)
1890 write (iout,*) 'too many processors for parallel prng'
1891 write (*,*) 'too many processors for parallel prng'
1899 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)