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 c if (modecalc.eq.8) then
32 c inquire (file="fort.40",exist=file_exist)
33 c 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
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 c include 'COMMON.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
73 c include 'COMMON.MAP'
74 include 'COMMON.HEADER'
75 c include 'COMMON.CSA'
76 include 'COMMON.CHAIN'
77 c include 'COMMON.MUCA'
79 include 'COMMON.FFIELD'
80 include 'COMMON.SETUP'
81 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
82 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
84 character*320 controlcard
89 read (INP,'(a)') titel
90 call card_concat(controlcard)
91 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
92 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
93 call reada(controlcard,'SEED',seed,0.0D0)
94 call random_init(seed)
95 C Set up the time limit (caution! The time must be input in minutes!)
96 read_cart=index(controlcard,'READ_CART').gt.0
97 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
98 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
99 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
100 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
101 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
102 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
103 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
104 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
105 call reada(controlcard,'DRMS',drms,0.1D0)
106 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
107 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
108 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
109 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
110 write (iout,'(a,f10.1)')'DRMS = ',drms
111 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
112 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
114 call readi(controlcard,'NZ_START',nz_start,0)
115 call readi(controlcard,'NZ_END',nz_end,0)
116 call readi(controlcard,'IZ_SC',iz_sc,0)
118 safety = 60.0d0*safety
121 call reada(controlcard,"T_BATH",t_bath,300.0d0)
122 minim=(index(controlcard,'MINIMIZE').gt.0)
123 dccart=(index(controlcard,'CART').gt.0)
124 overlapsc=(index(controlcard,'OVERLAP').gt.0)
125 overlapsc=.not.overlapsc
126 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
127 searchsc=.not.searchsc
128 sideadd=(index(controlcard,'SIDEADD').gt.0)
129 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
130 outpdb=(index(controlcard,'PDBOUT').gt.0)
131 outmol2=(index(controlcard,'MOL2OUT').gt.0)
132 pdbref=(index(controlcard,'PDBREF').gt.0)
133 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
134 indpdb=index(controlcard,'PDBSTART')
135 extconf=(index(controlcard,'EXTCONF').gt.0)
136 call readi(controlcard,'IPRINT',iprint,0)
137 call readi(controlcard,'MAXGEN',maxgen,10000)
138 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
139 call readi(controlcard,"KDIAG",kdiag,0)
140 call readi(controlcard,"RESCALE_MODE",rescale_mode,1)
141 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
142 & write (iout,*) "RESCALE_MODE",rescale_mode
143 split_ene=index(controlcard,'SPLIT_ENE').gt.0
144 if (index(controlcard,'REGULAR').gt.0.0D0) then
145 call reada(controlcard,'WEIDIS',weidis,0.1D0)
149 if (index(controlcard,'CHECKGRAD').gt.0) then
151 if (index(controlcard,'CART').gt.0) then
153 elseif (index(controlcard,'CARINT').gt.0) then
158 elseif (index(controlcard,'THREAD').gt.0) then
160 call readi(controlcard,'THREAD',nthread,0)
161 if (nthread.gt.0) then
162 call reada(controlcard,'WEIDIS',weidis,0.1D0)
165 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
166 stop 'Error termination in Read_Control.'
168 else if (index(controlcard,'MCMA').gt.0) then
170 else if (index(controlcard,'MCEE').gt.0) then
172 else if (index(controlcard,'MULTCONF').gt.0) then
174 else if (index(controlcard,'MAP').gt.0) then
176 call readi(controlcard,'MAP',nmap,0)
177 else if (index(controlcard,'CSA').gt.0) then
179 crc else if (index(controlcard,'ZSCORE').gt.0) then
181 crc ZSCORE is rm from UNRES, modecalc=9 is available
184 cfcm else if (index(controlcard,'MCMF').gt.0) then
186 else if (index(controlcard,'SOFTREG').gt.0) then
188 else if (index(controlcard,'CHECK_BOND').gt.0) then
190 else if (index(controlcard,'TEST').gt.0) then
192 else if (index(controlcard,'MD').gt.0) then
194 else if (index(controlcard,'RE ').gt.0) then
198 lmuca=index(controlcard,'MUCA').gt.0
199 call readi(controlcard,'MUCADYN',mucadyn,0)
200 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
201 if (lmuca .and. (me.eq.king .or. .not.out1file ))
203 write (iout,*) 'MUCADYN=',mucadyn
204 write (iout,*) 'MUCASMOOTH=',muca_smooth
207 iscode=index(controlcard,'ONE_LETTER')
208 indphi=index(controlcard,'PHI')
209 indback=index(controlcard,'BACK')
210 iranconf=index(controlcard,'RAND_CONF')
211 i2ndstr=index(controlcard,'USE_SEC_PRED')
212 gradout=index(controlcard,'GRADOUT').gt.0
213 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
215 if(me.eq.king.or..not.out1file)
216 & write (iout,'(2a)') diagmeth(kdiag),
217 & ' routine used to diagonalize matrices.'
220 c------------------------------------------------------------------------------
223 C Read molecular data.
225 implicit real*8 (a-h,o-z)
231 include 'COMMON.IOUNITS'
234 include 'COMMON.INTERACT'
235 include 'COMMON.LOCAL'
236 include 'COMMON.NAMES'
237 include 'COMMON.CHAIN'
238 include 'COMMON.FFIELD'
239 include 'COMMON.SBRIDGE'
240 include 'COMMON.HEADER'
241 include 'COMMON.CONTROL'
242 c include 'COMMON.DBASE'
243 c include 'COMMON.THREAD'
244 include 'COMMON.CONTACTS'
245 include 'COMMON.TORCNSTR'
246 include 'COMMON.TIME1'
247 include 'COMMON.BOUNDS'
248 c include 'COMMON.MD'
249 c include 'COMMON.REMD'
250 include 'COMMON.SETUP'
251 character*4 sequence(maxres)
253 double precision x(maxvar)
254 character*256 pdbfile
255 character*320 weightcard
256 character*80 weightcard_t,ucase
257 dimension itype_pdb(maxres)
258 common /pizda/ itype_pdb
259 logical seq_comp,fail
260 double precision energia(0:n_ene)
266 C Read weights of the subsequent energy terms.
268 call card_concat(weightcard)
269 call reada(weightcard,'WLONG',wlong,1.0D0)
270 call reada(weightcard,'WSC',wsc,wlong)
271 call reada(weightcard,'WSCP',wscp,wlong)
272 call reada(weightcard,'WELEC',welec,1.0D0)
273 call reada(weightcard,'WVDWPP',wvdwpp,welec)
274 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
275 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
276 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
277 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
278 call reada(weightcard,'WTURN3',wturn3,1.0D0)
279 call reada(weightcard,'WTURN4',wturn4,1.0D0)
280 call reada(weightcard,'WTURN6',wturn6,1.0D0)
281 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
282 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
283 call reada(weightcard,'WBOND',wbond,1.0D0)
284 call reada(weightcard,'WTOR',wtor,1.0D0)
285 call reada(weightcard,'WTORD',wtor_d,1.0D0)
286 call reada(weightcard,'WANG',wang,1.0D0)
287 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
288 call reada(weightcard,'SCAL14',scal14,0.4D0)
289 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
290 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
291 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
292 call reada(weightcard,'TEMP0',temp0,300.0d0)
293 if (index(weightcard,'SOFT').gt.0) ipot=6
294 C 12/1/95 Added weight for the multi-body term WCORR
295 call reada(weightcard,'WCORRH',wcorr,1.0D0)
296 if (wcorr4.gt.0.0d0) wcorr=wcorr4
317 if(me.eq.king.or..not.out1file)
318 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
319 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
321 10 format (/'Energy-term weights (unscaled):'//
322 & 'WSCC= ',f10.6,' (SC-SC)'/
323 & 'WSCP= ',f10.6,' (SC-p)'/
324 & 'WELEC= ',f10.6,' (p-p electr)'/
325 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
326 & 'WBOND= ',f10.6,' (stretching)'/
327 & 'WANG= ',f10.6,' (bending)'/
328 & 'WSCLOC= ',f10.6,' (SC local)'/
329 & 'WTOR= ',f10.6,' (torsional)'/
330 & 'WTORD= ',f10.6,' (double torsional)'/
331 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
332 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
333 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
334 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
335 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
336 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
337 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
338 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
339 & 'WTURN6= ',f10.6,' (turns, 6th order)')
340 if(me.eq.king.or..not.out1file)then
341 if (wcorr4.gt.0.0d0) then
342 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
343 & 'between contact pairs of peptide groups'
344 write (iout,'(2(a,f5.3/))')
345 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
346 & 'Range of quenching the correlation terms:',2*delt_corr
347 else if (wcorr.gt.0.0d0) then
348 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
349 & 'between contact pairs of peptide groups'
351 write (iout,'(a,f8.3)')
352 & 'Scaling factor of 1,4 SC-p interactions:',scal14
353 write (iout,'(a,f8.3)')
354 & 'General scaling factor of SC-p interactions:',scalscp
356 r0_corr=cutoff_corr-delt_corr
358 aad(i,1)=scalscp*aad(i,1)
359 aad(i,2)=scalscp*aad(i,2)
360 bad(i,1)=scalscp*bad(i,1)
361 bad(i,2)=scalscp*bad(i,2)
363 call rescale_weights(t_bath)
364 if(me.eq.king.or..not.out1file)
365 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
366 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
368 22 format (/'Energy-term weights (scaled):'//
369 & 'WSCC= ',f10.6,' (SC-SC)'/
370 & 'WSCP= ',f10.6,' (SC-p)'/
371 & 'WELEC= ',f10.6,' (p-p electr)'/
372 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
373 & 'WBOND= ',f10.6,' (stretching)'/
374 & 'WANG= ',f10.6,' (bending)'/
375 & 'WSCLOC= ',f10.6,' (SC local)'/
376 & 'WTOR= ',f10.6,' (torsional)'/
377 & 'WTORD= ',f10.6,' (double torsional)'/
378 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
379 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
380 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
381 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
382 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
383 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
384 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
385 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
386 & 'WTURN6= ',f10.6,' (turns, 6th order)')
387 if(me.eq.king.or..not.out1file)
388 & write (iout,*) "Reference temperature for weights calculation:",
390 call reada(weightcard,"D0CM",d0cm,3.78d0)
391 call reada(weightcard,"AKCM",akcm,15.1d0)
392 call reada(weightcard,"AKTH",akth,11.0d0)
393 call reada(weightcard,"AKCT",akct,12.0d0)
394 call reada(weightcard,"V1SS",v1ss,-1.08d0)
395 call reada(weightcard,"V2SS",v2ss,7.61d0)
396 call reada(weightcard,"V3SS",v3ss,13.7d0)
397 call reada(weightcard,"EBR",ebr,-5.50D0)
398 if(me.eq.king.or..not.out1file) then
399 write (iout,*) "Parameters of the SS-bond potential:"
400 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
402 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
403 write (iout,*) "EBR",ebr
404 print *,'indpdb=',indpdb,' pdbref=',pdbref
406 if (indpdb.gt.0 .or. pdbref) then
407 read(inp,'(a)') pdbfile
408 if(me.eq.king.or..not.out1file)
409 & write (iout,'(2a)') 'PDB data will be read from file ',
410 & pdbfile(:ilen(pdbfile))
411 open(ipdbin,file=pdbfile,status='old',err=33)
413 33 write (iout,'(a)') 'Error opening PDB file.'
416 c print *,'Begin reading pdb data'
418 c print *,'Finished reading pdb data'
419 if(me.eq.king.or..not.out1file)
420 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
421 & ' nstart_sup=',nstart_sup
423 itype_pdb(i)=itype(i)
427 nct=nstart_sup+nsup-1
428 c call contact(.false.,ncont_ref,icont_ref,co)
431 if(me.eq.king.or..not.out1file)
432 & write(iout,*)'Adding sidechains'
439 do while (fail.and.nsi.le.maxsi)
440 c call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
443 if(fail) write(iout,*)'Adding sidechain failed for res ',
444 & i,' after ',nsi,' trials'
449 if (indpdb.eq.0) then
450 C Read sequence if not taken from the pdb file.
452 c print *,'nres=',nres
453 if (iscode.gt.0) then
454 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
456 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
458 C Convert sequence to numeric code
460 itype(i)=rescode(i,sequence(i),iscode)
462 C Assign initial virtual bond lengths
468 vbld(i+nres)=dsc(itype(i))
469 vbld_inv(i+nres)=dsc_inv(itype(i))
470 c write (iout,*) "i",i," itype",itype(i),
471 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
475 c print '(20i4)',(itype(i),i=1,nres)
478 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
480 if (itype(i).eq.21) then
484 else if (itype(i+1).ne.20) then
486 else if (itype(i).ne.20) then
493 if(me.eq.king.or..not.out1file)then
494 write (iout,*) "ITEL"
496 write (iout,*) i,itype(i),itel(i)
498 print *,'Call Read_Bridge.'
501 C 8/13/98 Set limits to generating the dihedral angles
506 read (inp,*) ndih_constr
507 if (ndih_constr.gt.0) then
509 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
510 if(me.eq.king.or..not.out1file)then
512 & 'There are',ndih_constr,' constraints on phi angles.'
514 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
518 phi0(i)=deg2rad*phi0(i)
519 drange(i)=deg2rad*drange(i)
521 if(me.eq.king.or..not.out1file)
522 & write (iout,*) 'FTORS',ftors
525 phibound(1,ii) = phi0(i)-drange(i)
526 phibound(2,ii) = phi0(i)+drange(i)
533 write (iout,'(a)') 'Boundaries in phi angle sampling:'
535 write (iout,'(a3,i5,2f10.1)')
536 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
542 cd print *,'NNT=',NNT,' NCT=',NCT
543 if (itype(1).eq.21) nnt=2
544 if (itype(nres).eq.21) nct=nct-1
546 if(me.eq.king.or..not.out1file)
547 & write (iout,'(a,i3)') 'nsup=',nsup
549 if (nsup.le.(nct-nnt+1)) then
550 do i=0,nct-nnt+1-nsup
551 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
557 & 'Error - sequences to be superposed do not match.'
560 do i=0,nsup-(nct-nnt+1)
561 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
563 nstart_sup=nstart_sup+i
569 & 'Error - sequences to be superposed do not match.'
572 if (nsup.eq.0) nsup=nct-nnt
573 if (nstart_sup.eq.0) nstart_sup=nnt
574 if (nstart_seq.eq.0) nstart_seq=nnt
575 if(me.eq.king.or..not.out1file)
576 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
577 & ' nstart_seq=',nstart_seq
579 c--- Zscore rms -------
580 if (nz_start.eq.0) nz_start=nnt
581 if (nz_end.eq.0 .and. nsup.gt.0) then
583 else if (nz_end.eq.0) then
586 if(me.eq.king.or..not.out1file)then
587 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
588 write (iout,*) 'IZ_SC=',iz_sc
590 c----------------------
593 if (.not.pdbref) then
594 call read_angles(inp,*38)
596 38 write (iout,'(a)') 'Error reading reference structure.'
598 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
599 stop 'Error reading reference structure'
603 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
612 c call contact(.true.,ncont_ref,icont_ref,co)
614 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
616 if (constr_dist.gt.0) call read_dist_constr
617 c write (iout,*) "After read_dist_constr nhpb",nhpb
619 if(me.eq.king.or..not.out1file)
620 & write (iout,*) 'Contact order:',co
622 if(me.eq.king.or..not.out1file)
623 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
626 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
628 if(me.eq.king.or..not.out1file)
629 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
630 & icont_ref(1,i),' ',
631 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
635 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
636 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
637 & modecalc.ne.10) then
638 C If input structure hasn't been supplied from the PDB file read or generate
640 if (iranconf.eq.0 .and. .not. extconf) then
641 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
642 & write (iout,'(a)') 'Initial geometry will be read in.'
644 read(inp,'(8f10.5)',end=36,err=36)
645 & ((c(l,k),l=1,3),k=1,nres),
646 & ((c(l,k+nres),l=1,3),k=nnt,nct)
647 call int_from_cart1(.false.)
650 dc(j,i)=c(j,i+1)-c(j,i)
651 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
655 if (itype(i).ne.10) then
657 dc(j,i+nres)=c(j,i+nres)-c(j,i)
658 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
664 call read_angles(inp,*36)
667 36 write (iout,'(a)') 'Error reading angle file.'
669 call mpi_finalize( MPI_COMM_WORLD,IERR )
671 stop 'Error reading angle file.'
673 else if (extconf) then
674 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
675 & write (iout,'(a)') 'Extended chain initial geometry.'
677 theta(i)=90d0*deg2rad
683 alph(i)=110d0*deg2rad
686 omeg(i)=-120d0*deg2rad
689 if(me.eq.king.or..not.out1file)
690 & write (iout,'(a)') 'Random-generated initial geometry.'
694 if (me.eq.king .or. fg_rank.eq.0 .and. (
695 & modecalc.eq.12 .or. modecalc.eq.14) ) then
699 c call gen_rand_conf(itmp,*30)
701 30 write (iout,*) 'Failed to generate random conformation',
703 write (*,*) 'Processor:',me,
704 & ' Failed to generate random conformation',
714 write (iout,'(a,i3,a)') 'Processor:',me,
715 & ' error in generating random conformation.'
716 write (*,'(a,i3,a)') 'Processor:',me,
717 & ' error in generating random conformation.'
720 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
726 c call gen_rand_conf(itmp,*31)
728 31 write (iout,*) 'Failed to generate random conformation',
730 write (*,*) 'Failed to generate random conformation',
733 write (iout,'(a,i3,a)') 'Processor:',me,
734 & ' error in generating random conformation.'
735 write (*,'(a,i3,a)') 'Processor:',me,
736 & ' error in generating random conformation.'
741 elseif (modecalc.eq.4) then
742 read (inp,'(a)') intinname
743 open (intin,file=intinname,status='old',err=333)
744 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
745 & write (iout,'(a)') 'intinname',intinname
746 write (*,'(a)') 'Processor',myrank,' intinname',intinname
748 333 write (iout,'(2a)') 'Error opening angle file ',intinname
750 call MPI_Finalize(MPI_COMM_WORLD,IERR)
752 stop 'Error opening angle file.'
756 C Generate distance constraints, if the PDB structure is to be regularized.
757 if (nthread.gt.0) then
761 if (me.eq.king .or. .not. out1file)
763 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
764 write (iout,'(/a,i3,a)')
765 & 'The chain contains',ns,' disulfide-bridging cysteines.'
766 write (iout,'(20i4)') (iss(i),i=1,ns)
767 write (iout,'(/a/)') 'Pre-formed links are:'
773 if (me.eq.king.or..not.out1file)
774 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
775 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
780 c if (i2ndstr.gt.0) call secstrp2dihc
781 c call geom_to_var(nvar,x)
782 c call etotal(energia(0))
783 c call enerprint(energia(0))
784 c call briefout(0,etot)
786 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
787 cd write (iout,'(a)') 'Variable list:'
788 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
790 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
791 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
792 & 'Processor',myrank,': end reading molecular data.'
796 c--------------------------------------------------------------------------
797 logical function seq_comp(itypea,itypeb,length)
799 integer length,itypea(length),itypeb(length)
802 if (itypea(i).ne.itypeb(i)) then
810 c-----------------------------------------------------------------------------
811 subroutine read_bridge
812 C Read information about disulfide bridges.
813 implicit real*8 (a-h,o-z)
818 include 'COMMON.IOUNITS'
821 include 'COMMON.INTERACT'
822 include 'COMMON.LOCAL'
823 include 'COMMON.NAMES'
824 include 'COMMON.CHAIN'
825 include 'COMMON.FFIELD'
826 include 'COMMON.SBRIDGE'
827 include 'COMMON.HEADER'
828 include 'COMMON.CONTROL'
829 c include 'COMMON.DBASE'
830 c include 'COMMON.THREAD'
831 include 'COMMON.TIME1'
832 include 'COMMON.SETUP'
833 C Read bridging residues.
834 read (inp,*) ns,(iss(i),i=1,ns)
836 if(me.eq.king.or..not.out1file)
837 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
838 C Check whether the specified bridging residues are cystines.
840 if (itype(iss(i)).ne.1) then
841 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
842 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
843 & ' can form a disulfide bridge?!!!'
844 write (*,'(2a,i3,a)')
845 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
846 & ' can form a disulfide bridge?!!!'
848 call MPI_Finalize(MPI_COMM_WORLD,ierror)
853 C Read preformed bridges.
855 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
856 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
859 C Check if the residues involved in bridges are in the specified list of
863 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
864 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
865 write (iout,'(a,i3,a)') 'Disulfide pair',i,
866 & ' contains residues present in other pairs.'
867 write (*,'(a,i3,a)') 'Disulfide pair',i,
868 & ' contains residues present in other pairs.'
870 call MPI_Finalize(MPI_COMM_WORLD,ierror)
876 if (ihpb(i).eq.iss(j)) goto 10
878 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
881 if (jhpb(i).eq.iss(j)) goto 20
883 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
896 c----------------------------------------------------------------------------
897 subroutine read_x(kanal,*)
898 implicit real*8 (a-h,o-z)
902 include 'COMMON.CHAIN'
903 include 'COMMON.IOUNITS'
904 include 'COMMON.CONTROL'
905 include 'COMMON.LOCAL'
906 include 'COMMON.INTERACT'
907 c Read coordinates from input
909 read(kanal,'(8f10.5)',end=10,err=10)
910 & ((c(l,k),l=1,3),k=1,nres),
911 & ((c(l,k+nres),l=1,3),k=nnt,nct)
914 c(j,2*nres)=c(j,nres)
916 call int_from_cart1(.false.)
919 dc(j,i)=c(j,i+1)-c(j,i)
920 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
924 if (itype(i).ne.10) then
926 dc(j,i+nres)=c(j,i+nres)-c(j,i)
927 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
935 c------------------------------------------------------------------------------
937 implicit real*8 (a-h,o-z)
939 include 'COMMON.IOUNITS'
942 include 'COMMON.INTERACT'
943 include 'COMMON.LOCAL'
944 include 'COMMON.NAMES'
945 include 'COMMON.CHAIN'
946 include 'COMMON.FFIELD'
947 include 'COMMON.SBRIDGE'
948 include 'COMMON.HEADER'
949 include 'COMMON.CONTROL'
950 c include 'COMMON.DBASE'
951 c include 'COMMON.THREAD'
952 include 'COMMON.TIME1'
953 C Set up variable list.
959 if (itype(i).ne.10) then
961 ialph(i,1)=nvar+nside
965 if (indphi.gt.0) then
967 else if (indback.gt.0) then
972 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
975 c----------------------------------------------------------------------------
976 subroutine gen_dist_constr
977 C Generate CA distance constraints.
978 implicit real*8 (a-h,o-z)
980 include 'COMMON.IOUNITS'
983 include 'COMMON.INTERACT'
984 include 'COMMON.LOCAL'
985 include 'COMMON.NAMES'
986 include 'COMMON.CHAIN'
987 include 'COMMON.FFIELD'
988 include 'COMMON.SBRIDGE'
989 include 'COMMON.HEADER'
990 include 'COMMON.CONTROL'
991 c include 'COMMON.DBASE'
992 c include 'COMMON.THREAD'
993 include 'COMMON.TIME1'
994 dimension itype_pdb(maxres)
995 common /pizda/ itype_pdb
997 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
998 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
999 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1001 do i=nstart_sup,nstart_sup+nsup-1
1002 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1003 cd & ' seq_pdb', restyp(itype_pdb(i))
1004 do j=i+2,nstart_sup+nsup-1
1006 ihpb(nhpb)=i+nstart_seq-nstart_sup
1007 jhpb(nhpb)=j+nstart_seq-nstart_sup
1009 dhpb(nhpb)=dist(i,j)
1012 cd write (iout,'(a)') 'Distance constraints:'
1017 cd if (ii.gt.nres) then
1022 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1023 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1024 cd & dhpb(i),forcon(i)
1029 subroutine read_minim
1030 implicit real*8 (a-h,o-z)
1031 include 'DIMENSIONS'
1032 include 'COMMON.MINIM'
1033 include 'COMMON.IOUNITS'
1035 character*320 minimcard
1036 call card_concat(minimcard)
1037 call readi(minimcard,'MAXMIN',maxmin,2000)
1038 call readi(minimcard,'MAXFUN',maxfun,5000)
1039 call readi(minimcard,'MINMIN',minmin,maxmin)
1040 call readi(minimcard,'MINFUN',minfun,maxmin)
1041 call reada(minimcard,'TOLF',tolf,1.0D-2)
1042 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1043 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1044 & 'Options in energy minimization:'
1045 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1046 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1047 & 'MinMin:',MinMin,' MinFun:',MinFun,
1048 & ' TolF:',TolF,' RTolF:',RTolF
1051 c----------------------------------------------------------------------------
1052 subroutine read_angles(kanal,*)
1053 implicit real*8 (a-h,o-z)
1054 include 'DIMENSIONS'
1055 include 'COMMON.GEO'
1056 include 'COMMON.VAR'
1057 include 'COMMON.CHAIN'
1058 include 'COMMON.IOUNITS'
1059 include 'COMMON.CONTROL'
1060 c Read angles from input
1062 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1063 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1064 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1065 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1068 c 9/7/01 avoid 180 deg valence angle
1069 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1071 theta(i)=deg2rad*theta(i)
1072 phi(i)=deg2rad*phi(i)
1073 alph(i)=deg2rad*alph(i)
1074 omeg(i)=deg2rad*omeg(i)
1079 c----------------------------------------------------------------------------
1080 subroutine reada(rekord,lancuch,wartosc,default)
1082 character*(*) rekord,lancuch
1083 double precision wartosc,default
1086 iread=index(rekord,lancuch)
1087 if (iread.eq.0) then
1091 iread=iread+ilen(lancuch)+1
1092 read (rekord(iread:),*,err=10,end=10) wartosc
1097 c----------------------------------------------------------------------------
1098 subroutine readi(rekord,lancuch,wartosc,default)
1100 character*(*) rekord,lancuch
1101 integer wartosc,default
1104 iread=index(rekord,lancuch)
1105 if (iread.eq.0) then
1109 iread=iread+ilen(lancuch)+1
1110 read (rekord(iread:),*,err=10,end=10) wartosc
1115 c----------------------------------------------------------------------------
1116 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1119 integer tablica(dim),default
1120 character*(*) rekord,lancuch
1127 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1128 if (iread.eq.0) return
1129 iread=iread+ilen(lancuch)+1
1130 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1133 c----------------------------------------------------------------------------
1134 subroutine multreada(rekord,lancuch,tablica,dim,default)
1137 double precision tablica(dim),default
1138 character*(*) rekord,lancuch
1145 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1146 if (iread.eq.0) return
1147 iread=iread+ilen(lancuch)+1
1148 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1151 c----------------------------------------------------------------------------
1152 subroutine openunits
1153 implicit real*8 (a-h,o-z)
1154 include 'DIMENSIONS'
1157 character*16 form,nodename
1160 include 'COMMON.SETUP'
1161 include 'COMMON.IOUNITS'
1162 c include 'COMMON.MD'
1163 include 'COMMON.CONTROL'
1164 integer lenpre,lenpot,ilen,lentmp
1166 character*3 out1file_text,ucase
1169 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1170 call getenv_loc("PREFIX",prefix)
1172 call getenv_loc("POT",pot)
1173 call getenv_loc("DIRTMP",tmpdir)
1174 call getenv_loc("CURDIR",curdir)
1175 call getenv_loc("OUT1FILE",out1file_text)
1176 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1177 out1file_text=ucase(out1file_text)
1178 if (out1file_text(1:1).eq."Y") then
1181 out1file=fg_rank.gt.0
1186 if (lentmp.gt.0) then
1187 write (*,'(80(1h!))')
1188 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1189 write (*,'(80(1h!))')
1190 write (*,*)"All output files will be on node /tmp directory."
1192 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1193 if (me.eq.king) then
1194 write (*,*) "The master node is ",nodename
1195 else if (fg_rank.eq.0) then
1196 write (*,*) "I am the CG slave node ",nodename
1198 write (*,*) "I am the FG slave node ",nodename
1201 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1202 lenpre = lentmp+lenpre+1
1204 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1205 C Get the names and open the input files
1206 #if defined(WINIFL) || defined(WINPGI)
1207 open(1,file=pref_orig(:ilen(pref_orig))//
1208 & '.inp',status='old',readonly,shared)
1209 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1210 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1211 C Get parameter filenames and open the parameter files.
1212 call getenv_loc('BONDPAR',bondname)
1213 open (ibond,file=bondname,status='old',readonly,shared)
1214 call getenv_loc('THETPAR',thetname)
1215 open (ithep,file=thetname,status='old',readonly,shared)
1217 call getenv_loc('THETPARPDB',thetname_pdb)
1218 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1220 call getenv_loc('ROTPAR',rotname)
1221 open (irotam,file=rotname,status='old',readonly,shared)
1223 call getenv_loc('ROTPARPDB',rotname_pdb)
1224 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1226 call getenv_loc('TORPAR',torname)
1227 open (itorp,file=torname,status='old',readonly,shared)
1228 call getenv_loc('TORDPAR',tordname)
1229 open (itordp,file=tordname,status='old',readonly,shared)
1230 call getenv_loc('FOURIER',fouriername)
1231 open (ifourier,file=fouriername,status='old',readonly,shared)
1232 call getenv_loc('ELEPAR',elename)
1233 open (ielep,file=elename,status='old',readonly,shared)
1234 call getenv_loc('SIDEPAR',sidename)
1235 open (isidep,file=sidename,status='old',readonly,shared)
1236 #elif (defined CRAY) || (defined AIX)
1237 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1239 c print *,"Processor",myrank," opened file 1"
1240 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1241 c print *,"Processor",myrank," opened file 9"
1242 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1243 C Get parameter filenames and open the parameter files.
1244 call getenv_loc('BONDPAR',bondname)
1245 open (ibond,file=bondname,status='old',action='read')
1246 c print *,"Processor",myrank," opened file IBOND"
1247 call getenv_loc('THETPAR',thetname)
1248 open (ithep,file=thetname,status='old',action='read')
1249 c print *,"Processor",myrank," opened file ITHEP"
1251 call getenv_loc('THETPARPDB',thetname_pdb)
1252 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1254 call getenv_loc('ROTPAR',rotname)
1255 open (irotam,file=rotname,status='old',action='read')
1256 c print *,"Processor",myrank," opened file IROTAM"
1258 call getenv_loc('ROTPARPDB',rotname_pdb)
1259 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1261 call getenv_loc('TORPAR',torname)
1262 open (itorp,file=torname,status='old',action='read')
1263 c print *,"Processor",myrank," opened file ITORP"
1264 call getenv_loc('TORDPAR',tordname)
1265 open (itordp,file=tordname,status='old',action='read')
1266 c print *,"Processor",myrank," opened file ITORDP"
1267 call getenv_loc('SCCORPAR',sccorname)
1268 open (isccor,file=sccorname,status='old',action='read')
1269 c print *,"Processor",myrank," opened file ISCCOR"
1270 call getenv_loc('FOURIER',fouriername)
1271 open (ifourier,file=fouriername,status='old',action='read')
1272 c print *,"Processor",myrank," opened file IFOURIER"
1273 call getenv_loc('ELEPAR',elename)
1274 open (ielep,file=elename,status='old',action='read')
1275 c print *,"Processor",myrank," opened file IELEP"
1276 call getenv_loc('SIDEPAR',sidename)
1277 open (isidep,file=sidename,status='old',action='read')
1278 c print *,"Processor",myrank," opened file ISIDEP"
1279 c print *,"Processor",myrank," opened parameter files"
1281 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1282 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1283 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1284 C Get parameter filenames and open the parameter files.
1285 call getenv_loc('BONDPAR',bondname)
1286 open (ibond,file=bondname,status='old')
1287 call getenv_loc('THETPAR',thetname)
1288 open (ithep,file=thetname,status='old')
1290 call getenv_loc('THETPARPDB',thetname_pdb)
1291 open (ithep_pdb,file=thetname_pdb,status='old')
1293 call getenv_loc('ROTPAR',rotname)
1294 open (irotam,file=rotname,status='old')
1296 call getenv_loc('ROTPARPDB',rotname_pdb)
1297 open (irotam_pdb,file=rotname_pdb,status='old')
1299 call getenv_loc('TORPAR',torname)
1300 open (itorp,file=torname,status='old')
1301 call getenv_loc('TORDPAR',tordname)
1302 open (itordp,file=tordname,status='old')
1303 call getenv_loc('SCCORPAR',sccorname)
1304 open (isccor,file=sccorname,status='old')
1305 call getenv_loc('FOURIER',fouriername)
1306 open (ifourier,file=fouriername,status='old')
1307 call getenv_loc('ELEPAR',elename)
1308 open (ielep,file=elename,status='old')
1309 call getenv_loc('SIDEPAR',sidename)
1310 open (isidep,file=sidename,status='old')
1312 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1314 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1315 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1316 C Get parameter filenames and open the parameter files.
1317 call getenv_loc('BONDPAR',bondname)
1318 open (ibond,file=bondname,status='old',readonly)
1319 call getenv_loc('THETPAR',thetname)
1320 open (ithep,file=thetname,status='old',readonly)
1322 call getenv_loc('THETPARPDB',thetname_pdb)
1323 print *,"thetname_pdb ",thetname_pdb
1324 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1325 print *,ithep_pdb," opened"
1327 call getenv_loc('ROTPAR',rotname)
1328 open (irotam,file=rotname,status='old',readonly)
1330 call getenv_loc('ROTPARPDB',rotname_pdb)
1331 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1333 call getenv_loc('TORPAR',torname)
1334 open (itorp,file=torname,status='old',readonly)
1335 call getenv_loc('TORDPAR',tordname)
1336 open (itordp,file=tordname,status='old',readonly)
1337 call getenv_loc('SCCORPAR',sccorname)
1338 open (isccor,file=sccorname,status='old',readonly)
1339 call getenv_loc('FOURIER',fouriername)
1340 open (ifourier,file=fouriername,status='old',readonly)
1341 call getenv_loc('ELEPAR',elename)
1342 open (ielep,file=elename,status='old',readonly)
1343 call getenv_loc('SIDEPAR',sidename)
1344 open (isidep,file=sidename,status='old',readonly)
1348 C 8/9/01 In the newest version SCp interaction constants are read from a file
1349 C Use -DOLDSCP to use hard-coded constants instead.
1351 call getenv_loc('SCPPAR',scpname)
1352 #if defined(WINIFL) || defined(WINPGI)
1353 open (iscpp,file=scpname,status='old',readonly,shared)
1354 #elif (defined CRAY) || (defined AIX)
1355 open (iscpp,file=scpname,status='old',action='read')
1357 open (iscpp,file=scpname,status='old')
1359 open (iscpp,file=scpname,status='old',readonly)
1362 call getenv_loc('PATTERN',patname)
1363 #if defined(WINIFL) || defined(WINPGI)
1364 open (icbase,file=patname,status='old',readonly,shared)
1365 #elif (defined CRAY) || (defined AIX)
1366 open (icbase,file=patname,status='old',action='read')
1368 open (icbase,file=patname,status='old')
1370 open (icbase,file=patname,status='old',readonly)
1373 C Open output file only for CG processes
1374 c print *,"Processor",myrank," fg_rank",fg_rank
1375 if (fg_rank.eq.0) then
1377 if (nodes.eq.1) then
1380 npos = dlog10(dfloat(nodes-1))+1
1382 if (npos.lt.3) npos=3
1383 write (liczba,'(i1)') npos
1384 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1386 write (liczba,form) me
1387 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1388 & liczba(:ilen(liczba))
1389 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1391 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1393 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1394 & liczba(:ilen(liczba))//'.mol2'
1395 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1396 & liczba(:ilen(liczba))//'.stat'
1398 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1399 & //liczba(:ilen(liczba))//'.stat')
1400 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1403 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1404 c & liczba(:ilen(liczba))//'.const'
1409 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1410 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1411 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1412 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1413 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1415 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1417 rest2name=prefix(:ilen(prefix))//'.rst'
1419 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1422 #if defined(AIX) || defined(PGI)
1423 if (me.eq.king .or. .not. out1file)
1424 & open(iout,file=outname,status='unknown')
1426 if (fg_rank.gt.0) then
1427 write (liczba,'(i3.3)') myrank/nfgtasks
1428 write (ll,'(bz,i3.3)') fg_rank
1429 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1434 open(igeom,file=intname,status='unknown',position='append')
1435 open(ipdb,file=pdbname,status='unknown')
1436 open(imol2,file=mol2name,status='unknown')
1437 open(istat,file=statname,status='unknown',position='append')
1439 c1out open(iout,file=outname,status='unknown')
1442 if (me.eq.king .or. .not.out1file)
1443 & open(iout,file=outname,status='unknown')
1445 if (fg_rank.gt.0) then
1446 write (liczba,'(i3.3)') myrank/nfgtasks
1447 write (ll,'(bz,i3.3)') fg_rank
1448 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1453 open(igeom,file=intname,status='unknown',access='append')
1454 open(ipdb,file=pdbname,status='unknown')
1455 open(imol2,file=mol2name,status='unknown')
1456 open(istat,file=statname,status='unknown',access='append')
1458 c1out open(iout,file=outname,status='unknown')
1461 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1462 csa_seed=prefix(:lenpre)//'.CSA.seed'
1463 csa_history=prefix(:lenpre)//'.CSA.history'
1464 csa_bank=prefix(:lenpre)//'.CSA.bank'
1465 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1466 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1467 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1468 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1469 csa_int=prefix(:lenpre)//'.int'
1470 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1471 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1472 csa_in=prefix(:lenpre)//'.CSA.in'
1473 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1476 write (iout,'(80(1h-))')
1477 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1478 write (iout,'(80(1h-))')
1479 write (iout,*) "Input file : ",
1480 & pref_orig(:ilen(pref_orig))//'.inp'
1481 write (iout,*) "Output file : ",
1482 & outname(:ilen(outname))
1484 write (iout,*) "Sidechain potential file : ",
1485 & sidename(:ilen(sidename))
1487 write (iout,*) "SCp potential file : ",
1488 & scpname(:ilen(scpname))
1490 write (iout,*) "Electrostatic potential file : ",
1491 & elename(:ilen(elename))
1492 write (iout,*) "Cumulant coefficient file : ",
1493 & fouriername(:ilen(fouriername))
1494 write (iout,*) "Torsional parameter file : ",
1495 & torname(:ilen(torname))
1496 write (iout,*) "Double torsional parameter file : ",
1497 & tordname(:ilen(tordname))
1498 write (iout,*) "SCCOR parameter file : ",
1499 & sccorname(:ilen(sccorname))
1500 write (iout,*) "Bond & inertia constant file : ",
1501 & bondname(:ilen(bondname))
1502 write (iout,*) "Bending parameter file : ",
1503 & thetname(:ilen(thetname))
1504 write (iout,*) "Rotamer parameter file : ",
1505 & rotname(:ilen(rotname))
1506 write (iout,*) "Threading database : ",
1507 & patname(:ilen(patname))
1509 &write (iout,*)" DIRTMP : ",
1511 write (iout,'(80(1h-))')
1515 c----------------------------------------------------------------------------
1516 subroutine card_concat(card)
1517 implicit real*8 (a-h,o-z)
1518 include 'DIMENSIONS'
1519 include 'COMMON.IOUNITS'
1521 character*80 karta,ucase
1523 read (inp,'(a)') karta
1526 do while (karta(80:80).eq.'&')
1527 card=card(:ilen(card)+1)//karta(:79)
1528 read (inp,'(a)') karta
1531 card=card(:ilen(card)+1)//karta
1535 subroutine read_dist_constr
1536 implicit real*8 (a-h,o-z)
1537 include 'DIMENSIONS'
1541 include 'COMMON.SETUP'
1542 include 'COMMON.CONTROL'
1543 include 'COMMON.CHAIN'
1544 include 'COMMON.IOUNITS'
1545 include 'COMMON.SBRIDGE'
1546 integer ifrag_(2,100),ipair_(2,100)
1547 double precision wfrag_(100),wpair_(100)
1548 character*500 controlcard
1549 c write (iout,*) "Calling read_dist_constr"
1550 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1552 call card_concat(controlcard)
1553 call readi(controlcard,"NFRAG",nfrag_,0)
1554 call readi(controlcard,"NPAIR",npair_,0)
1555 call readi(controlcard,"NDIST",ndist_,0)
1556 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1557 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1558 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1559 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1560 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1561 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1562 c write (iout,*) "IFRAG"
1564 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1566 c write (iout,*) "IPAIR"
1568 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1572 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1573 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1574 & ifrag_(2,i)=nstart_sup+nsup-1
1575 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1577 if (wfrag_(i).gt.0.0d0) then
1578 do j=ifrag_(1,i),ifrag_(2,i)-1
1579 do k=j+1,ifrag_(2,i)
1580 write (iout,*) "j",j," k",k
1582 if (constr_dist.eq.1) then
1587 forcon(nhpb)=wfrag_(i)
1588 else if (constr_dist.eq.2) then
1589 if (ddjk.le.dist_cut) then
1594 forcon(nhpb)=wfrag_(i)
1601 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1604 if (.not.out1file .or. me.eq.king)
1605 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1606 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1608 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1609 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1616 if (wpair_(i).gt.0.0d0) then
1624 do j=ifrag_(1,ii),ifrag_(2,ii)
1625 do k=ifrag_(1,jj),ifrag_(2,jj)
1629 forcon(nhpb)=wpair_(i)
1630 dhpb(nhpb)=dist(j,k)
1632 if (.not.out1file .or. me.eq.king)
1633 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1634 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1636 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1637 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1644 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1645 if (forcon(nhpb+1).gt.0.0d0) then
1647 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1649 if (.not.out1file .or. me.eq.king)
1650 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1651 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1653 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1654 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1661 c-------------------------------------------------------------------------------
1663 subroutine flush(iu)
1668 subroutine flush(iu)
1673 c------------------------------------------------------------------------------
1674 subroutine copy_to_tmp(source)
1675 include "DIMENSIONS"
1676 include "COMMON.IOUNITS"
1677 character*(*) source
1678 character* 256 tmpfile
1682 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1683 inquire(file=tmpfile,exist=ex)
1685 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1686 & " to temporary directory..."
1687 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1688 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1692 c------------------------------------------------------------------------------
1693 subroutine move_from_tmp(source)
1694 include "DIMENSIONS"
1695 include "COMMON.IOUNITS"
1696 character*(*) source
1699 write (*,*) "Moving ",source(:ilen(source)),
1700 & " from temporary directory to working directory"
1701 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1702 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1705 c------------------------------------------------------------------------------
1706 subroutine random_init(seed)
1708 C Initialize random number generator
1710 implicit real*8 (a-h,o-z)
1711 include 'DIMENSIONS'
1717 logical OKRandom, prng_restart
1719 integer iseed_array(4)
1721 include 'COMMON.IOUNITS'
1722 include 'COMMON.TIME1'
1723 c include 'COMMON.THREAD'
1724 include 'COMMON.SBRIDGE'
1725 include 'COMMON.CONTROL'
1726 include 'COMMON.MCM'
1727 c include 'COMMON.MAP'
1728 include 'COMMON.HEADER'
1729 c include 'COMMON.CSA'
1730 include 'COMMON.CHAIN'
1731 c include 'COMMON.MUCA'
1732 c include 'COMMON.MD'
1733 include 'COMMON.FFIELD'
1734 include 'COMMON.SETUP'
1735 iseed=-dint(dabs(seed))
1736 if (iseed.eq.0) then
1737 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1738 & 'Random seed undefined. The program will stop.'
1739 write (*,'(/80(1h*)/20x,a/80(1h*))')
1740 & 'Random seed undefined. The program will stop.'
1742 call mpi_finalize(mpi_comm_world,ierr)
1744 stop 'Bad random seed.'
1747 if (fg_rank.eq.0) then
1751 if(me.eq.king .or. .not. out1file)
1752 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1753 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1754 OKRandom = prng_restart(me,iseedi8)
1757 tmp=65536.0d0**(4-i)
1758 iseed_array(i) = dint(seed/tmp)
1759 seed=seed-iseed_array(i)*tmp
1761 if(me.eq.king .or. .not. out1file)
1762 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1763 & (iseed_array(i),i=1,4)
1764 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1765 & (iseed_array(i),i=1,4)
1766 OKRandom = prng_restart(me,iseed_array)
1769 c r1=ran_number(0.0D0,1.0D0)
1770 if(me.eq.king .or. .not. out1file)
1771 & write (iout,*) 'ran_num',r1
1772 if (r1.lt.0.0d0) OKRandom=.false.
1774 if (.not.OKRandom) then
1775 write (iout,*) 'PRNG IS NOT WORKING!!!'
1776 print *,'PRNG IS NOT WORKING!!!'
1779 call mpi_abort(mpi_comm_world,error_msg,ierr)
1782 write (iout,*) 'too many processors for parallel prng'
1783 write (*,*) 'too many processors for parallel prng'
1791 c write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)