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 include 'COMMON.SBRIDGE'
70 include 'COMMON.CONTROL'
72 include 'COMMON.HEADER'
73 include 'COMMON.CHAIN'
74 include 'COMMON.FFIELD'
75 include 'COMMON.SETUP'
76 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
77 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
79 character*320 controlcard
84 read (INP,'(a)') titel
85 call card_concat(controlcard)
86 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
87 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
88 call reada(controlcard,'SEED',seed,0.0D0)
89 call random_init(seed)
90 C Set up the time limit (caution! The time must be input in minutes!)
91 read_cart=index(controlcard,'READ_CART').gt.0
92 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
93 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
94 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
95 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
96 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
97 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
98 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
99 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
100 call reada(controlcard,'DRMS',drms,0.1D0)
101 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
102 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
103 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
104 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
105 write (iout,'(a,f10.1)')'DRMS = ',drms
106 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
107 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
109 call readi(controlcard,'NZ_START',nz_start,0)
110 call readi(controlcard,'NZ_END',nz_end,0)
111 call readi(controlcard,'IZ_SC',iz_sc,0)
113 safety = 60.0d0*safety
116 call reada(controlcard,"T_BATH",t_bath,300.0d0)
117 minim=(index(controlcard,'MINIMIZE').gt.0)
118 dccart=(index(controlcard,'CART').gt.0)
119 overlapsc=(index(controlcard,'OVERLAP').gt.0)
120 overlapsc=.not.overlapsc
121 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
122 searchsc=.not.searchsc
123 sideadd=(index(controlcard,'SIDEADD').gt.0)
124 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
125 outpdb=(index(controlcard,'PDBOUT').gt.0)
126 outmol2=(index(controlcard,'MOL2OUT').gt.0)
127 pdbref=(index(controlcard,'PDBREF').gt.0)
128 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
129 indpdb=index(controlcard,'PDBSTART')
130 extconf=(index(controlcard,'EXTCONF').gt.0)
131 call readi(controlcard,'IPRINT',iprint,0)
132 call readi(controlcard,'MAXGEN',maxgen,10000)
133 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
134 call readi(controlcard,"KDIAG",kdiag,0)
135 call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
136 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
137 & write (iout,*) "RESCALE_MODE",rescale_mode
138 split_ene=index(controlcard,'SPLIT_ENE').gt.0
139 if (index(controlcard,'REGULAR').gt.0.0D0) then
140 call reada(controlcard,'WEIDIS',weidis,0.1D0)
144 if (index(controlcard,'CHECKGRAD').gt.0) then
146 if (index(controlcard,'CART').gt.0) then
148 elseif (index(controlcard,'CARINT').gt.0) then
153 elseif (index(controlcard,'THREAD').gt.0) then
155 call readi(controlcard,'THREAD',nthread,0)
156 if (nthread.gt.0) then
157 call reada(controlcard,'WEIDIS',weidis,0.1D0)
160 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
161 stop 'Error termination in Read_Control.'
163 else if (index(controlcard,'MCMA').gt.0) then
165 else if (index(controlcard,'MCEE').gt.0) then
167 else if (index(controlcard,'MULTCONF').gt.0) then
169 else if (index(controlcard,'MAP').gt.0) then
171 call readi(controlcard,'MAP',nmap,0)
172 else if (index(controlcard,'CSA').gt.0) then
174 crc else if (index(controlcard,'ZSCORE').gt.0) then
176 crc ZSCORE is rm from UNRES, modecalc=9 is available
179 cfcm else if (index(controlcard,'MCMF').gt.0) then
181 else if (index(controlcard,'SOFTREG').gt.0) then
183 else if (index(controlcard,'CHECK_BOND').gt.0) then
185 else if (index(controlcard,'TEST').gt.0) then
187 else if (index(controlcard,'MD').gt.0) then
189 else if (index(controlcard,'RE ').gt.0) then
193 lmuca=index(controlcard,'MUCA').gt.0
194 call readi(controlcard,'MUCADYN',mucadyn,0)
195 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
196 if (lmuca .and. (me.eq.king .or. .not.out1file ))
198 write (iout,*) 'MUCADYN=',mucadyn
199 write (iout,*) 'MUCASMOOTH=',muca_smooth
202 iscode=index(controlcard,'ONE_LETTER')
203 indphi=index(controlcard,'PHI')
204 indback=index(controlcard,'BACK')
205 iranconf=index(controlcard,'RAND_CONF')
206 i2ndstr=index(controlcard,'USE_SEC_PRED')
207 gradout=index(controlcard,'GRADOUT').gt.0
208 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
210 if(me.eq.king.or..not.out1file)
211 & write (iout,'(2a)') diagmeth(kdiag),
212 & ' routine used to diagonalize matrices.'
215 c------------------------------------------------------------------------------
218 C Read molecular data.
220 implicit real*8 (a-h,o-z)
226 include 'COMMON.IOUNITS'
229 include 'COMMON.INTERACT'
230 include 'COMMON.LOCAL'
231 include 'COMMON.NAMES'
232 include 'COMMON.CHAIN'
233 include 'COMMON.FFIELD'
234 include 'COMMON.SBRIDGE'
235 include 'COMMON.HEADER'
236 include 'COMMON.CONTROL'
237 c include 'COMMON.DBASE'
238 c include 'COMMON.THREAD'
239 include 'COMMON.CONTACTS'
240 include 'COMMON.TORCNSTR'
241 include 'COMMON.TIME1'
242 include 'COMMON.BOUNDS'
243 c include 'COMMON.MD'
244 c include 'COMMON.REMD'
245 include 'COMMON.SETUP'
246 character*4 sequence(maxres)
248 double precision x(maxvar)
249 character*256 pdbfile
250 character*320 weightcard
251 character*80 weightcard_t,ucase
252 dimension itype_pdb(maxres)
253 common /pizda/ itype_pdb
254 logical seq_comp,fail
255 double precision energia(0:n_ene)
261 C Read weights of the subsequent energy terms.
263 call card_concat(weightcard)
264 call reada(weightcard,'WLONG',wlong,1.0D0)
265 call reada(weightcard,'WSC',wsc,wlong)
266 call reada(weightcard,'WSCP',wscp,wlong)
267 call reada(weightcard,'WELEC',welec,1.0D0)
268 call reada(weightcard,'WVDWPP',wvdwpp,welec)
269 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
270 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
271 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
272 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
273 call reada(weightcard,'WTURN3',wturn3,1.0D0)
274 call reada(weightcard,'WTURN4',wturn4,1.0D0)
275 call reada(weightcard,'WTURN6',wturn6,1.0D0)
276 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
277 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
278 call reada(weightcard,'WBOND',wbond,1.0D0)
279 call reada(weightcard,'WTOR',wtor,1.0D0)
280 call reada(weightcard,'WTORD',wtor_d,1.0D0)
281 call reada(weightcard,'WANG',wang,1.0D0)
282 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
283 call reada(weightcard,'SCAL14',scal14,0.4D0)
284 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
285 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
286 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
287 call reada(weightcard,'TEMP0',temp0,300.0d0)
288 if (index(weightcard,'SOFT').gt.0) ipot=6
289 C 12/1/95 Added weight for the multi-body term WCORR
290 call reada(weightcard,'WCORRH',wcorr,1.0D0)
291 if (wcorr4.gt.0.0d0) wcorr=wcorr4
312 if(me.eq.king.or..not.out1file)
313 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
314 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
316 10 format (/'Energy-term weights (unscaled):'//
317 & 'WSCC= ',f10.6,' (SC-SC)'/
318 & 'WSCP= ',f10.6,' (SC-p)'/
319 & 'WELEC= ',f10.6,' (p-p electr)'/
320 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
321 & 'WBOND= ',f10.6,' (stretching)'/
322 & 'WANG= ',f10.6,' (bending)'/
323 & 'WSCLOC= ',f10.6,' (SC local)'/
324 & 'WTOR= ',f10.6,' (torsional)'/
325 & 'WTORD= ',f10.6,' (double torsional)'/
326 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
327 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
328 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
329 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
330 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
331 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
332 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
333 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
334 & 'WTURN6= ',f10.6,' (turns, 6th order)')
335 if(me.eq.king.or..not.out1file)then
336 if (wcorr4.gt.0.0d0) then
337 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
338 & 'between contact pairs of peptide groups'
339 write (iout,'(2(a,f5.3/))')
340 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
341 & 'Range of quenching the correlation terms:',2*delt_corr
342 else if (wcorr.gt.0.0d0) then
343 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
344 & 'between contact pairs of peptide groups'
346 write (iout,'(a,f8.3)')
347 & 'Scaling factor of 1,4 SC-p interactions:',scal14
348 write (iout,'(a,f8.3)')
349 & 'General scaling factor of SC-p interactions:',scalscp
351 r0_corr=cutoff_corr-delt_corr
353 aad(i,1)=scalscp*aad(i,1)
354 aad(i,2)=scalscp*aad(i,2)
355 bad(i,1)=scalscp*bad(i,1)
356 bad(i,2)=scalscp*bad(i,2)
358 c call rescale_weights(t_bath)
359 if(me.eq.king.or..not.out1file)
360 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
361 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
363 22 format (/'Energy-term weights (scaled):'//
364 & 'WSCC= ',f10.6,' (SC-SC)'/
365 & 'WSCP= ',f10.6,' (SC-p)'/
366 & 'WELEC= ',f10.6,' (p-p electr)'/
367 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
368 & 'WBOND= ',f10.6,' (stretching)'/
369 & 'WANG= ',f10.6,' (bending)'/
370 & 'WSCLOC= ',f10.6,' (SC local)'/
371 & 'WTOR= ',f10.6,' (torsional)'/
372 & 'WTORD= ',f10.6,' (double torsional)'/
373 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
374 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
375 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
376 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
377 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
378 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
379 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
380 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
381 & 'WTURN6= ',f10.6,' (turns, 6th order)')
382 if(me.eq.king.or..not.out1file)
383 & write (iout,*) "Reference temperature for weights calculation:",
385 call reada(weightcard,"D0CM",d0cm,3.78d0)
386 call reada(weightcard,"AKCM",akcm,15.1d0)
387 call reada(weightcard,"AKTH",akth,11.0d0)
388 call reada(weightcard,"AKCT",akct,12.0d0)
389 call reada(weightcard,"V1SS",v1ss,-1.08d0)
390 call reada(weightcard,"V2SS",v2ss,7.61d0)
391 call reada(weightcard,"V3SS",v3ss,13.7d0)
392 call reada(weightcard,"EBR",ebr,-5.50D0)
393 if(me.eq.king.or..not.out1file) then
394 write (iout,*) "Parameters of the SS-bond potential:"
395 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
397 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
398 write (iout,*) "EBR",ebr
399 print *,'indpdb=',indpdb,' pdbref=',pdbref
401 if (indpdb.gt.0 .or. pdbref) then
402 read(inp,'(a)') pdbfile
403 if(me.eq.king.or..not.out1file)
404 & write (iout,'(2a)') 'PDB data will be read from file ',
405 & pdbfile(:ilen(pdbfile))
406 open(ipdbin,file=pdbfile,status='old',err=33)
408 33 write (iout,'(a)') 'Error opening PDB file.'
411 c print *,'Begin reading pdb data'
413 c print *,'Finished reading pdb data'
414 if(me.eq.king.or..not.out1file)
415 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
416 & ' nstart_sup=',nstart_sup
418 itype_pdb(i)=itype(i)
422 nct=nstart_sup+nsup-1
423 c call contact(.false.,ncont_ref,icont_ref,co)
426 C Following 2 lines for diagnostics; comment out if not needed
427 write (iout,*) "Before sideadd"
429 if(me.eq.king.or..not.out1file)
430 & write(iout,*)'Adding sidechains'
437 do while (fail.and.nsi.le.maxsi)
438 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
441 if(fail) write(iout,*)'Adding sidechain failed for res ',
442 & i,' after ',nsi,' trials'
445 C 9/29/12 Adam: Recalculate coordinates with new side chain positions
448 C Following 2 lines for diagnostics; comment out if not needed
449 write (iout,*) "After sideadd"
452 if (indpdb.eq.0) then
453 C Read sequence if not taken from the pdb file.
455 c print *,'nres=',nres
456 if (iscode.gt.0) then
457 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
459 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
461 C Convert sequence to numeric code
463 itype(i)=rescode(i,sequence(i),iscode)
465 C Assign initial virtual bond lengths
471 vbld(i+nres)=dsc(itype(i))
472 vbld_inv(i+nres)=dsc_inv(itype(i))
473 c write (iout,*) "i",i," itype",itype(i),
474 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
478 c print '(20i4)',(itype(i),i=1,nres)
481 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
483 if (itype(i).eq.21) then
487 else if (itype(i+1).ne.20) then
489 else if (itype(i).ne.20) then
496 if(me.eq.king.or..not.out1file)then
497 write (iout,*) "ITEL"
499 write (iout,*) i,itype(i),itel(i)
501 print *,'Call Read_Bridge.'
504 C 8/13/98 Set limits to generating the dihedral angles
509 read (inp,*) ndih_constr
510 if (ndih_constr.gt.0) then
512 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
513 if(me.eq.king.or..not.out1file)then
515 & 'There are',ndih_constr,' constraints on phi angles.'
517 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
521 phi0(i)=deg2rad*phi0(i)
522 drange(i)=deg2rad*drange(i)
524 if(me.eq.king.or..not.out1file)
525 & write (iout,*) 'FTORS',ftors
528 phibound(1,ii) = phi0(i)-drange(i)
529 phibound(2,ii) = phi0(i)+drange(i)
536 write (iout,'(a)') 'Boundaries in phi angle sampling:'
538 write (iout,'(a3,i5,2f10.1)')
539 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
545 cd print *,'NNT=',NNT,' NCT=',NCT
546 if (itype(1).eq.21) nnt=2
547 if (itype(nres).eq.21) nct=nct-1
549 if(me.eq.king.or..not.out1file)
550 & write (iout,'(a,i3)') 'nsup=',nsup
552 if (nsup.le.(nct-nnt+1)) then
553 do i=0,nct-nnt+1-nsup
554 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
560 & 'Error - sequences to be superposed do not match.'
563 do i=0,nsup-(nct-nnt+1)
564 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
566 nstart_sup=nstart_sup+i
572 & 'Error - sequences to be superposed do not match.'
575 if (nsup.eq.0) nsup=nct-nnt
576 if (nstart_sup.eq.0) nstart_sup=nnt
577 if (nstart_seq.eq.0) nstart_seq=nnt
578 if(me.eq.king.or..not.out1file)
579 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
580 & ' nstart_seq=',nstart_seq
582 c--- Zscore rms -------
583 if (nz_start.eq.0) nz_start=nnt
584 if (nz_end.eq.0 .and. nsup.gt.0) then
586 else if (nz_end.eq.0) then
589 if(me.eq.king.or..not.out1file)then
590 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
591 write (iout,*) 'IZ_SC=',iz_sc
593 c----------------------
596 if (.not.pdbref) then
597 call read_angles(inp,*38)
599 38 write (iout,'(a)') 'Error reading reference structure.'
601 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
602 stop 'Error reading reference structure'
606 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
615 c call contact(.true.,ncont_ref,icont_ref,co)
617 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
619 if (constr_dist.gt.0) call read_dist_constr
620 c write (iout,*) "After read_dist_constr nhpb",nhpb
622 if(me.eq.king.or..not.out1file)
623 & write (iout,*) 'Contact order:',co
625 if(me.eq.king.or..not.out1file)
626 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
629 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
631 if(me.eq.king.or..not.out1file)
632 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
633 & icont_ref(1,i),' ',
634 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
638 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
639 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
640 & modecalc.ne.10) then
641 C If input structure hasn't been supplied from the PDB file read or generate
643 if (iranconf.eq.0 .and. .not. extconf) then
644 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
645 & write (iout,'(a)') 'Initial geometry will be read in.'
647 read(inp,'(8f10.5)',end=36,err=36)
648 & ((c(l,k),l=1,3),k=1,nres),
649 & ((c(l,k+nres),l=1,3),k=nnt,nct)
650 call int_from_cart1(.false.)
653 dc(j,i)=c(j,i+1)-c(j,i)
654 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
658 if (itype(i).ne.10) then
660 dc(j,i+nres)=c(j,i+nres)-c(j,i)
661 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
667 call read_angles(inp,*36)
670 36 write (iout,'(a)') 'Error reading angle file.'
672 call mpi_finalize( MPI_COMM_WORLD,IERR )
674 stop 'Error reading angle file.'
676 else if (extconf) then
677 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
678 & write (iout,'(a)') 'Extended chain initial geometry.'
680 theta(i)=90d0*deg2rad
686 alph(i)=110d0*deg2rad
689 omeg(i)=-120d0*deg2rad
692 if(me.eq.king.or..not.out1file)
693 & write (iout,'(a)') 'Random-generated initial geometry.'
697 if (me.eq.king .or. fg_rank.eq.0 .and. (
698 & modecalc.eq.12 .or. modecalc.eq.14) ) then
702 call gen_rand_conf(itmp,*30)
704 30 write (iout,*) 'Failed to generate random conformation',
706 write (*,*) 'Processor:',me,
707 & ' Failed to generate random conformation',
717 write (iout,'(a,i3,a)') 'Processor:',me,
718 & ' error in generating random conformation.'
719 write (*,'(a,i3,a)') 'Processor:',me,
720 & ' error in generating random conformation.'
723 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
729 call gen_rand_conf(itmp,*31)
731 31 write (iout,*) 'Failed to generate random conformation',
733 write (*,*) 'Failed to generate random conformation',
736 write (iout,'(a,i3,a)') 'Processor:',me,
737 & ' error in generating random conformation.'
738 write (*,'(a,i3,a)') 'Processor:',me,
739 & ' error in generating random conformation.'
744 elseif (modecalc.eq.4) then
745 read (inp,'(a)') intinname
746 open (intin,file=intinname,status='old',err=333)
747 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
748 & write (iout,'(a)') 'intinname',intinname
749 write (*,'(a)') 'Processor',myrank,' intinname',intinname
751 333 write (iout,'(2a)') 'Error opening angle file ',intinname
753 call MPI_Finalize(MPI_COMM_WORLD,IERR)
755 stop 'Error opening angle file.'
759 C Generate distance constraints, if the PDB structure is to be regularized.
760 c if (nthread.gt.0) then
761 c call read_threadbase
764 if (me.eq.king .or. .not. out1file)
766 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
767 write (iout,'(/a,i3,a)')
768 & 'The chain contains',ns,' disulfide-bridging cysteines.'
769 write (iout,'(20i4)') (iss(i),i=1,ns)
770 write (iout,'(/a/)') 'Pre-formed links are:'
776 if (me.eq.king.or..not.out1file)
777 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
778 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
783 c if (i2ndstr.gt.0) call secstrp2dihc
784 c call geom_to_var(nvar,x)
785 c call etotal(energia(0))
786 c call enerprint(energia(0))
787 c call briefout(0,etot)
789 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
790 cd write (iout,'(a)') 'Variable list:'
791 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
793 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
794 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
795 & 'Processor',myrank,': end reading molecular data.'
799 c--------------------------------------------------------------------------
800 logical function seq_comp(itypea,itypeb,length)
802 integer length,itypea(length),itypeb(length)
805 if (itypea(i).ne.itypeb(i)) then
813 c-----------------------------------------------------------------------------
814 subroutine read_bridge
815 C Read information about disulfide bridges.
816 implicit real*8 (a-h,o-z)
821 include 'COMMON.IOUNITS'
824 include 'COMMON.INTERACT'
825 include 'COMMON.LOCAL'
826 include 'COMMON.NAMES'
827 include 'COMMON.CHAIN'
828 include 'COMMON.FFIELD'
829 include 'COMMON.SBRIDGE'
830 include 'COMMON.HEADER'
831 include 'COMMON.CONTROL'
832 c include 'COMMON.DBASE'
833 c include 'COMMON.THREAD'
834 include 'COMMON.TIME1'
835 include 'COMMON.SETUP'
836 C Read bridging residues.
837 read (inp,*) ns,(iss(i),i=1,ns)
839 if(me.eq.king.or..not.out1file)
840 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
841 C Check whether the specified bridging residues are cystines.
843 if (itype(iss(i)).ne.1) then
844 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
845 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
846 & ' can form a disulfide bridge?!!!'
847 write (*,'(2a,i3,a)')
848 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
849 & ' can form a disulfide bridge?!!!'
851 call MPI_Finalize(MPI_COMM_WORLD,ierror)
856 C Read preformed bridges.
858 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
859 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
862 C Check if the residues involved in bridges are in the specified list of
866 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
867 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
868 write (iout,'(a,i3,a)') 'Disulfide pair',i,
869 & ' contains residues present in other pairs.'
870 write (*,'(a,i3,a)') 'Disulfide pair',i,
871 & ' contains residues present in other pairs.'
873 call MPI_Finalize(MPI_COMM_WORLD,ierror)
879 if (ihpb(i).eq.iss(j)) goto 10
881 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
884 if (jhpb(i).eq.iss(j)) goto 20
886 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
899 c----------------------------------------------------------------------------
900 subroutine read_x(kanal,*)
901 implicit real*8 (a-h,o-z)
905 include 'COMMON.CHAIN'
906 include 'COMMON.IOUNITS'
907 include 'COMMON.CONTROL'
908 include 'COMMON.LOCAL'
909 include 'COMMON.INTERACT'
910 c Read coordinates from input
912 read(kanal,'(8f10.5)',end=10,err=10)
913 & ((c(l,k),l=1,3),k=1,nres),
914 & ((c(l,k+nres),l=1,3),k=nnt,nct)
917 c(j,2*nres)=c(j,nres)
919 call int_from_cart1(.false.)
922 dc(j,i)=c(j,i+1)-c(j,i)
923 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
927 if (itype(i).ne.10) then
929 dc(j,i+nres)=c(j,i+nres)-c(j,i)
930 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
938 c------------------------------------------------------------------------------
940 implicit real*8 (a-h,o-z)
942 include 'COMMON.IOUNITS'
945 include 'COMMON.INTERACT'
946 include 'COMMON.LOCAL'
947 include 'COMMON.NAMES'
948 include 'COMMON.CHAIN'
949 include 'COMMON.FFIELD'
950 include 'COMMON.SBRIDGE'
951 include 'COMMON.HEADER'
952 include 'COMMON.CONTROL'
953 c include 'COMMON.DBASE'
954 c include 'COMMON.THREAD'
955 include 'COMMON.TIME1'
956 C Set up variable list.
962 if (itype(i).ne.10) then
964 ialph(i,1)=nvar+nside
968 if (indphi.gt.0) then
970 else if (indback.gt.0) then
975 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
978 c----------------------------------------------------------------------------
979 subroutine gen_dist_constr
980 C Generate CA distance constraints.
981 implicit real*8 (a-h,o-z)
983 include 'COMMON.IOUNITS'
986 include 'COMMON.INTERACT'
987 include 'COMMON.LOCAL'
988 include 'COMMON.NAMES'
989 include 'COMMON.CHAIN'
990 include 'COMMON.FFIELD'
991 include 'COMMON.SBRIDGE'
992 include 'COMMON.HEADER'
993 include 'COMMON.CONTROL'
994 c include 'COMMON.DBASE'
995 c include 'COMMON.THREAD'
996 include 'COMMON.TIME1'
997 dimension itype_pdb(maxres)
998 common /pizda/ itype_pdb
1000 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1001 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1002 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1004 do i=nstart_sup,nstart_sup+nsup-1
1005 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1006 cd & ' seq_pdb', restyp(itype_pdb(i))
1007 do j=i+2,nstart_sup+nsup-1
1009 ihpb(nhpb)=i+nstart_seq-nstart_sup
1010 jhpb(nhpb)=j+nstart_seq-nstart_sup
1012 dhpb(nhpb)=dist(i,j)
1015 cd write (iout,'(a)') 'Distance constraints:'
1020 cd if (ii.gt.nres) then
1025 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1026 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1027 cd & dhpb(i),forcon(i)
1032 subroutine read_minim
1033 implicit real*8 (a-h,o-z)
1034 include 'DIMENSIONS'
1035 include 'COMMON.MINIM'
1036 include 'COMMON.IOUNITS'
1038 character*320 minimcard
1039 call card_concat(minimcard)
1040 call readi(minimcard,'MAXMIN',maxmin,2000)
1041 call readi(minimcard,'MAXFUN',maxfun,5000)
1042 call readi(minimcard,'MINMIN',minmin,maxmin)
1043 call readi(minimcard,'MINFUN',minfun,maxmin)
1044 call reada(minimcard,'TOLF',tolf,1.0D-2)
1045 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1046 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1047 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1048 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1049 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1050 & 'Options in energy minimization:'
1051 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1052 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1053 & 'MinMin:',MinMin,' MinFun:',MinFun,
1054 & ' TolF:',TolF,' RTolF:',RTolF
1057 c----------------------------------------------------------------------------
1058 subroutine read_angles(kanal,*)
1059 implicit real*8 (a-h,o-z)
1060 include 'DIMENSIONS'
1061 include 'COMMON.GEO'
1062 include 'COMMON.VAR'
1063 include 'COMMON.CHAIN'
1064 include 'COMMON.IOUNITS'
1065 include 'COMMON.CONTROL'
1066 c Read angles from input
1068 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1069 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1070 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1071 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1074 c 9/7/01 avoid 180 deg valence angle
1075 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1077 theta(i)=deg2rad*theta(i)
1078 phi(i)=deg2rad*phi(i)
1079 alph(i)=deg2rad*alph(i)
1080 omeg(i)=deg2rad*omeg(i)
1085 c----------------------------------------------------------------------------
1086 subroutine reada(rekord,lancuch,wartosc,default)
1088 character*(*) rekord,lancuch
1089 double precision wartosc,default
1092 iread=index(rekord,lancuch)
1093 if (iread.eq.0) then
1097 iread=iread+ilen(lancuch)+1
1098 read (rekord(iread:),*,err=10,end=10) wartosc
1103 c----------------------------------------------------------------------------
1104 subroutine readi(rekord,lancuch,wartosc,default)
1106 character*(*) rekord,lancuch
1107 integer wartosc,default
1110 iread=index(rekord,lancuch)
1111 if (iread.eq.0) then
1115 iread=iread+ilen(lancuch)+1
1116 read (rekord(iread:),*,err=10,end=10) wartosc
1121 c----------------------------------------------------------------------------
1122 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1125 integer tablica(dim),default
1126 character*(*) rekord,lancuch
1133 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1134 if (iread.eq.0) return
1135 iread=iread+ilen(lancuch)+1
1136 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1139 c----------------------------------------------------------------------------
1140 subroutine multreada(rekord,lancuch,tablica,dim,default)
1143 double precision tablica(dim),default
1144 character*(*) rekord,lancuch
1151 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1152 if (iread.eq.0) return
1153 iread=iread+ilen(lancuch)+1
1154 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1157 c----------------------------------------------------------------------------
1158 subroutine openunits
1159 implicit real*8 (a-h,o-z)
1160 include 'DIMENSIONS'
1163 character*16 form,nodename
1166 include 'COMMON.SETUP'
1167 include 'COMMON.IOUNITS'
1168 c include 'COMMON.MD'
1169 include 'COMMON.CONTROL'
1170 integer lenpre,lenpot,ilen,lentmp
1172 character*3 out1file_text,ucase
1175 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1176 call getenv_loc("PREFIX",prefix)
1178 call getenv_loc("POT",pot)
1179 call getenv_loc("DIRTMP",tmpdir)
1180 call getenv_loc("CURDIR",curdir)
1181 call getenv_loc("OUT1FILE",out1file_text)
1182 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1183 out1file_text=ucase(out1file_text)
1184 if (out1file_text(1:1).eq."Y") then
1187 out1file=fg_rank.gt.0
1192 if (lentmp.gt.0) then
1193 write (*,'(80(1h!))')
1194 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1195 write (*,'(80(1h!))')
1196 write (*,*)"All output files will be on node /tmp directory."
1198 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1199 if (me.eq.king) then
1200 write (*,*) "The master node is ",nodename
1201 else if (fg_rank.eq.0) then
1202 write (*,*) "I am the CG slave node ",nodename
1204 write (*,*) "I am the FG slave node ",nodename
1207 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1208 lenpre = lentmp+lenpre+1
1210 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1211 C Get the names and open the input files
1212 #if defined(WINIFL) || defined(WINPGI)
1213 open(1,file=pref_orig(:ilen(pref_orig))//
1214 & '.inp',status='old',readonly,shared)
1215 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1216 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1217 C Get parameter filenames and open the parameter files.
1218 call getenv_loc('BONDPAR',bondname)
1219 open (ibond,file=bondname,status='old',readonly,shared)
1220 call getenv_loc('THETPAR',thetname)
1221 open (ithep,file=thetname,status='old',readonly,shared)
1223 call getenv_loc('THETPARPDB',thetname_pdb)
1224 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1226 call getenv_loc('ROTPAR',rotname)
1227 open (irotam,file=rotname,status='old',readonly,shared)
1229 call getenv_loc('ROTPARPDB',rotname_pdb)
1230 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1232 call getenv_loc('TORPAR',torname)
1233 open (itorp,file=torname,status='old',readonly,shared)
1234 call getenv_loc('TORDPAR',tordname)
1235 open (itordp,file=tordname,status='old',readonly,shared)
1236 call getenv_loc('FOURIER',fouriername)
1237 open (ifourier,file=fouriername,status='old',readonly,shared)
1238 call getenv_loc('ELEPAR',elename)
1239 open (ielep,file=elename,status='old',readonly,shared)
1240 call getenv_loc('SIDEPAR',sidename)
1241 open (isidep,file=sidename,status='old',readonly,shared)
1242 #elif (defined CRAY) || (defined AIX)
1243 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1245 c print *,"Processor",myrank," opened file 1"
1246 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1247 c print *,"Processor",myrank," opened file 9"
1248 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1249 C Get parameter filenames and open the parameter files.
1250 call getenv_loc('BONDPAR',bondname)
1251 open (ibond,file=bondname,status='old',action='read')
1252 c print *,"Processor",myrank," opened file IBOND"
1253 call getenv_loc('THETPAR',thetname)
1254 open (ithep,file=thetname,status='old',action='read')
1255 c print *,"Processor",myrank," opened file ITHEP"
1257 call getenv_loc('THETPARPDB',thetname_pdb)
1258 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1260 call getenv_loc('ROTPAR',rotname)
1261 open (irotam,file=rotname,status='old',action='read')
1262 c print *,"Processor",myrank," opened file IROTAM"
1264 call getenv_loc('ROTPARPDB',rotname_pdb)
1265 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1267 call getenv_loc('TORPAR',torname)
1268 open (itorp,file=torname,status='old',action='read')
1269 c print *,"Processor",myrank," opened file ITORP"
1270 call getenv_loc('TORDPAR',tordname)
1271 open (itordp,file=tordname,status='old',action='read')
1272 c print *,"Processor",myrank," opened file ITORDP"
1273 call getenv_loc('SCCORPAR',sccorname)
1274 open (isccor,file=sccorname,status='old',action='read')
1275 c print *,"Processor",myrank," opened file ISCCOR"
1276 call getenv_loc('FOURIER',fouriername)
1277 open (ifourier,file=fouriername,status='old',action='read')
1278 c print *,"Processor",myrank," opened file IFOURIER"
1279 call getenv_loc('ELEPAR',elename)
1280 open (ielep,file=elename,status='old',action='read')
1281 c print *,"Processor",myrank," opened file IELEP"
1282 call getenv_loc('SIDEPAR',sidename)
1283 open (isidep,file=sidename,status='old',action='read')
1284 c print *,"Processor",myrank," opened file ISIDEP"
1285 c print *,"Processor",myrank," opened parameter files"
1287 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1288 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1289 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1290 C Get parameter filenames and open the parameter files.
1291 call getenv_loc('BONDPAR',bondname)
1292 open (ibond,file=bondname,status='old')
1293 call getenv_loc('THETPAR',thetname)
1294 open (ithep,file=thetname,status='old')
1296 call getenv_loc('THETPARPDB',thetname_pdb)
1297 open (ithep_pdb,file=thetname_pdb,status='old')
1299 call getenv_loc('ROTPAR',rotname)
1300 open (irotam,file=rotname,status='old')
1302 call getenv_loc('ROTPARPDB',rotname_pdb)
1303 open (irotam_pdb,file=rotname_pdb,status='old')
1305 call getenv_loc('TORPAR',torname)
1306 open (itorp,file=torname,status='old')
1307 call getenv_loc('TORDPAR',tordname)
1308 open (itordp,file=tordname,status='old')
1309 call getenv_loc('SCCORPAR',sccorname)
1310 open (isccor,file=sccorname,status='old')
1311 call getenv_loc('FOURIER',fouriername)
1312 open (ifourier,file=fouriername,status='old')
1313 call getenv_loc('ELEPAR',elename)
1314 open (ielep,file=elename,status='old')
1315 call getenv_loc('SIDEPAR',sidename)
1316 open (isidep,file=sidename,status='old')
1318 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1320 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1321 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1322 C Get parameter filenames and open the parameter files.
1323 call getenv_loc('BONDPAR',bondname)
1324 open (ibond,file=bondname,status='old',readonly)
1325 call getenv_loc('THETPAR',thetname)
1326 open (ithep,file=thetname,status='old',readonly)
1328 call getenv_loc('THETPARPDB',thetname_pdb)
1329 print *,"thetname_pdb ",thetname_pdb
1330 open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1331 print *,ithep_pdb," opened"
1333 call getenv_loc('ROTPAR',rotname)
1334 open (irotam,file=rotname,status='old',readonly)
1336 call getenv_loc('ROTPARPDB',rotname_pdb)
1337 open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1339 call getenv_loc('TORPAR',torname)
1340 open (itorp,file=torname,status='old',readonly)
1341 call getenv_loc('TORDPAR',tordname)
1342 open (itordp,file=tordname,status='old',readonly)
1343 call getenv_loc('SCCORPAR',sccorname)
1344 open (isccor,file=sccorname,status='old',readonly)
1345 call getenv_loc('FOURIER',fouriername)
1346 open (ifourier,file=fouriername,status='old',readonly)
1347 call getenv_loc('ELEPAR',elename)
1348 open (ielep,file=elename,status='old',readonly)
1349 call getenv_loc('SIDEPAR',sidename)
1350 open (isidep,file=sidename,status='old',readonly)
1354 C 8/9/01 In the newest version SCp interaction constants are read from a file
1355 C Use -DOLDSCP to use hard-coded constants instead.
1357 call getenv_loc('SCPPAR',scpname)
1358 #if defined(WINIFL) || defined(WINPGI)
1359 open (iscpp,file=scpname,status='old',readonly,shared)
1360 #elif (defined CRAY) || (defined AIX)
1361 open (iscpp,file=scpname,status='old',action='read')
1363 open (iscpp,file=scpname,status='old')
1365 open (iscpp,file=scpname,status='old',readonly)
1368 call getenv_loc('PATTERN',patname)
1369 #if defined(WINIFL) || defined(WINPGI)
1370 open (icbase,file=patname,status='old',readonly,shared)
1371 #elif (defined CRAY) || (defined AIX)
1372 open (icbase,file=patname,status='old',action='read')
1374 open (icbase,file=patname,status='old')
1376 open (icbase,file=patname,status='old',readonly)
1379 C Open output file only for CG processes
1380 c print *,"Processor",myrank," fg_rank",fg_rank
1381 if (fg_rank.eq.0) then
1383 if (nodes.eq.1) then
1386 npos = dlog10(dfloat(nodes-1))+1
1388 if (npos.lt.3) npos=3
1389 write (liczba,'(i1)') npos
1390 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1392 write (liczba,form) me
1393 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1394 & liczba(:ilen(liczba))
1395 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1397 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1399 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1400 & liczba(:ilen(liczba))//'.mol2'
1401 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1402 & liczba(:ilen(liczba))//'.stat'
1404 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1405 & //liczba(:ilen(liczba))//'.stat')
1406 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1409 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1410 c & liczba(:ilen(liczba))//'.const'
1415 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1416 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1417 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1418 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1419 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1421 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1423 rest2name=prefix(:ilen(prefix))//'.rst'
1425 c qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1428 #if defined(AIX) || defined(PGI)
1429 if (me.eq.king .or. .not. out1file)
1430 & open(iout,file=outname,status='unknown')
1432 if (fg_rank.gt.0) then
1433 write (liczba,'(i3.3)') myrank/nfgtasks
1434 write (ll,'(bz,i3.3)') fg_rank
1435 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1440 open(igeom,file=intname,status='unknown',position='append')
1441 open(ipdb,file=pdbname,status='unknown')
1442 open(imol2,file=mol2name,status='unknown')
1443 open(istat,file=statname,status='unknown',position='append')
1445 c1out open(iout,file=outname,status='unknown')
1448 if (me.eq.king .or. .not.out1file)
1449 & open(iout,file=outname,status='unknown')
1451 if (fg_rank.gt.0) then
1452 write (liczba,'(i3.3)') myrank/nfgtasks
1453 write (ll,'(bz,i3.3)') fg_rank
1454 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1459 open(igeom,file=intname,status='unknown',access='append')
1460 open(ipdb,file=pdbname,status='unknown')
1461 open(imol2,file=mol2name,status='unknown')
1462 open(istat,file=statname,status='unknown',access='append')
1464 c1out open(iout,file=outname,status='unknown')
1467 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1468 csa_seed=prefix(:lenpre)//'.CSA.seed'
1469 csa_history=prefix(:lenpre)//'.CSA.history'
1470 csa_bank=prefix(:lenpre)//'.CSA.bank'
1471 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1472 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1473 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1474 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1475 csa_int=prefix(:lenpre)//'.int'
1476 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1477 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1478 csa_in=prefix(:lenpre)//'.CSA.in'
1479 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1482 write (iout,'(80(1h-))')
1483 write (iout,'(30x,a)') "FILE ASSIGNMENT"
1484 write (iout,'(80(1h-))')
1485 write (iout,*) "Input file : ",
1486 & pref_orig(:ilen(pref_orig))//'.inp'
1487 write (iout,*) "Output file : ",
1488 & outname(:ilen(outname))
1490 write (iout,*) "Sidechain potential file : ",
1491 & sidename(:ilen(sidename))
1493 write (iout,*) "SCp potential file : ",
1494 & scpname(:ilen(scpname))
1496 write (iout,*) "Electrostatic potential file : ",
1497 & elename(:ilen(elename))
1498 write (iout,*) "Cumulant coefficient file : ",
1499 & fouriername(:ilen(fouriername))
1500 write (iout,*) "Torsional parameter file : ",
1501 & torname(:ilen(torname))
1502 write (iout,*) "Double torsional parameter file : ",
1503 & tordname(:ilen(tordname))
1504 write (iout,*) "SCCOR parameter file : ",
1505 & sccorname(:ilen(sccorname))
1506 write (iout,*) "Bond & inertia constant file : ",
1507 & bondname(:ilen(bondname))
1508 write (iout,*) "Bending parameter file : ",
1509 & thetname(:ilen(thetname))
1510 write (iout,*) "Rotamer parameter file : ",
1511 & rotname(:ilen(rotname))
1512 write (iout,*) "Threading database : ",
1513 & patname(:ilen(patname))
1515 &write (iout,*)" DIRTMP : ",
1517 write (iout,'(80(1h-))')
1521 c----------------------------------------------------------------------------
1522 subroutine card_concat(card)
1523 implicit real*8 (a-h,o-z)
1524 include 'DIMENSIONS'
1525 include 'COMMON.IOUNITS'
1527 character*80 karta,ucase
1529 read (inp,'(a)') karta
1532 do while (karta(80:80).eq.'&')
1533 card=card(:ilen(card)+1)//karta(:79)
1534 read (inp,'(a)') karta
1537 card=card(:ilen(card)+1)//karta
1541 subroutine read_dist_constr
1542 implicit real*8 (a-h,o-z)
1543 include 'DIMENSIONS'
1547 include 'COMMON.SETUP'
1548 include 'COMMON.CONTROL'
1549 include 'COMMON.CHAIN'
1550 include 'COMMON.IOUNITS'
1551 include 'COMMON.SBRIDGE'
1552 integer ifrag_(2,100),ipair_(2,100)
1553 double precision wfrag_(100),wpair_(100)
1554 character*500 controlcard
1555 c write (iout,*) "Calling read_dist_constr"
1556 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1558 call card_concat(controlcard)
1559 call readi(controlcard,"NFRAG",nfrag_,0)
1560 call readi(controlcard,"NPAIR",npair_,0)
1561 call readi(controlcard,"NDIST",ndist_,0)
1562 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1563 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1564 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1565 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1566 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1567 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1568 c write (iout,*) "IFRAG"
1570 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1572 c write (iout,*) "IPAIR"
1574 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1578 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1579 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1580 & ifrag_(2,i)=nstart_sup+nsup-1
1581 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1583 if (wfrag_(i).gt.0.0d0) then
1584 do j=ifrag_(1,i),ifrag_(2,i)-1
1585 do k=j+1,ifrag_(2,i)
1586 write (iout,*) "j",j," k",k
1588 if (constr_dist.eq.1) then
1593 forcon(nhpb)=wfrag_(i)
1594 else if (constr_dist.eq.2) then
1595 if (ddjk.le.dist_cut) then
1600 forcon(nhpb)=wfrag_(i)
1607 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1610 if (.not.out1file .or. me.eq.king)
1611 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1612 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1614 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1615 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1622 if (wpair_(i).gt.0.0d0) then
1630 do j=ifrag_(1,ii),ifrag_(2,ii)
1631 do k=ifrag_(1,jj),ifrag_(2,jj)
1635 forcon(nhpb)=wpair_(i)
1636 dhpb(nhpb)=dist(j,k)
1638 if (.not.out1file .or. me.eq.king)
1639 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1640 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1642 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1643 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1650 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1651 if (forcon(nhpb+1).gt.0.0d0) then
1653 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1655 if (.not.out1file .or. me.eq.king)
1656 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1657 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1659 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1660 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1667 c-------------------------------------------------------------------------------
1669 subroutine flush(iu)
1674 subroutine flush(iu)
1679 c------------------------------------------------------------------------------
1680 subroutine copy_to_tmp(source)
1681 include "DIMENSIONS"
1682 include "COMMON.IOUNITS"
1683 character*(*) source
1684 character* 256 tmpfile
1688 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1689 inquire(file=tmpfile,exist=ex)
1691 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1692 & " to temporary directory..."
1693 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1694 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1698 c------------------------------------------------------------------------------
1699 subroutine move_from_tmp(source)
1700 include "DIMENSIONS"
1701 include "COMMON.IOUNITS"
1702 character*(*) source
1705 write (*,*) "Moving ",source(:ilen(source)),
1706 & " from temporary directory to working directory"
1707 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1708 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1711 c------------------------------------------------------------------------------
1712 subroutine random_init(seed)
1714 C Initialize random number generator
1716 implicit real*8 (a-h,o-z)
1717 include 'DIMENSIONS'
1723 logical OKRandom, prng_restart
1725 integer iseed_array(4)
1727 include 'COMMON.IOUNITS'
1728 include 'COMMON.TIME1'
1729 c include 'COMMON.THREAD'
1730 include 'COMMON.SBRIDGE'
1731 include 'COMMON.CONTROL'
1732 include 'COMMON.MCM'
1733 c include 'COMMON.MAP'
1734 include 'COMMON.HEADER'
1735 c include 'COMMON.CSA'
1736 include 'COMMON.CHAIN'
1737 c include 'COMMON.MUCA'
1738 c include 'COMMON.MD'
1739 include 'COMMON.FFIELD'
1740 include 'COMMON.SETUP'
1741 iseed=-dint(dabs(seed))
1742 if (iseed.eq.0) then
1743 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1744 & 'Random seed undefined. The program will stop.'
1745 write (*,'(/80(1h*)/20x,a/80(1h*))')
1746 & 'Random seed undefined. The program will stop.'
1748 call mpi_finalize(mpi_comm_world,ierr)
1750 stop 'Bad random seed.'
1753 if (fg_rank.eq.0) then
1757 if(me.eq.king .or. .not. out1file)
1758 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1759 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1760 OKRandom = prng_restart(me,iseedi8)
1763 tmp=65536.0d0**(4-i)
1764 iseed_array(i) = dint(seed/tmp)
1765 seed=seed-iseed_array(i)*tmp
1767 if(me.eq.king .or. .not. out1file)
1768 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1769 & (iseed_array(i),i=1,4)
1770 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1771 & (iseed_array(i),i=1,4)
1772 OKRandom = prng_restart(me,iseed_array)
1775 c r1=ran_number(0.0D0,1.0D0)
1776 if(me.eq.king .or. .not. out1file)
1777 & write (iout,*) 'ran_num',r1
1778 if (r1.lt.0.0d0) OKRandom=.false.
1780 if (.not.OKRandom) then
1781 write (iout,*) 'PRNG IS NOT WORKING!!!'
1782 print *,'PRNG IS NOT WORKING!!!'
1785 call mpi_abort(mpi_comm_world,error_msg,ierr)
1788 write (iout,*) 'too many processors for parallel prng'
1789 write (*,*) 'too many processors for parallel prng'
1797 c write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)