2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 call readi(controlcard,'SYM',symetr,1)
102 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
103 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
104 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
105 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
106 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
107 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
108 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
109 call reada(controlcard,'DRMS',drms,0.1D0)
110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
111 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
112 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
113 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
114 write (iout,'(a,f10.1)')'DRMS = ',drms
115 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
116 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
118 call readi(controlcard,'NZ_START',nz_start,0)
119 call readi(controlcard,'NZ_END',nz_end,0)
120 call readi(controlcard,'IZ_SC',iz_sc,0)
122 safety = 60.0d0*safety
125 call reada(controlcard,"T_BATH",t_bath,300.0d0)
126 minim=(index(controlcard,'MINIMIZE').gt.0)
127 dccart=(index(controlcard,'CART').gt.0)
128 overlapsc=(index(controlcard,'OVERLAP').gt.0)
129 overlapsc=.not.overlapsc
130 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
131 searchsc=.not.searchsc
132 sideadd=(index(controlcard,'SIDEADD').gt.0)
133 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
134 outpdb=(index(controlcard,'PDBOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 call readi(controlcard,'IPRINT',iprint,0)
141 call readi(controlcard,'MAXGEN',maxgen,10000)
142 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
143 call readi(controlcard,"KDIAG",kdiag,0)
144 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
145 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
146 & write (iout,*) "RESCALE_MODE",rescale_mode
147 split_ene=index(controlcard,'SPLIT_ENE').gt.0
148 if (index(controlcard,'REGULAR').gt.0.0D0) then
149 call reada(controlcard,'WEIDIS',weidis,0.1D0)
153 if (index(controlcard,'CHECKGRAD').gt.0) then
155 if (index(controlcard,'CART').gt.0) then
157 elseif (index(controlcard,'CARINT').gt.0) then
162 elseif (index(controlcard,'THREAD').gt.0) then
164 call readi(controlcard,'THREAD',nthread,0)
165 if (nthread.gt.0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
169 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
170 stop 'Error termination in Read_Control.'
172 else if (index(controlcard,'MCMA').gt.0) then
174 else if (index(controlcard,'MCEE').gt.0) then
176 else if (index(controlcard,'MULTCONF').gt.0) then
178 else if (index(controlcard,'MAP').gt.0) then
180 call readi(controlcard,'MAP',nmap,0)
181 else if (index(controlcard,'CSA').gt.0) then
183 crc else if (index(controlcard,'ZSCORE').gt.0) then
185 crc ZSCORE is rm from UNRES, modecalc=9 is available
188 cfcm else if (index(controlcard,'MCMF').gt.0) then
190 else if (index(controlcard,'SOFTREG').gt.0) then
192 else if (index(controlcard,'CHECK_BOND').gt.0) then
194 else if (index(controlcard,'TEST').gt.0) then
196 else if (index(controlcard,'MD').gt.0) then
198 else if (index(controlcard,'RE ').gt.0) then
202 lmuca=index(controlcard,'MUCA').gt.0
203 call readi(controlcard,'MUCADYN',mucadyn,0)
204 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
205 if (lmuca .and. (me.eq.king .or. .not.out1file ))
207 write (iout,*) 'MUCADYN=',mucadyn
208 write (iout,*) 'MUCASMOOTH=',muca_smooth
211 iscode=index(controlcard,'ONE_LETTER')
212 indphi=index(controlcard,'PHI')
213 indback=index(controlcard,'BACK')
214 iranconf=index(controlcard,'RAND_CONF')
215 i2ndstr=index(controlcard,'USE_SEC_PRED')
216 gradout=index(controlcard,'GRADOUT').gt.0
217 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
218 C DISTCHAINMAX become obsolete for periodic boundry condition
219 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
220 C Reading the dimensions of box in x,y,z coordinates
221 call reada(controlcard,'BOXX',boxxsize,100.0d0)
222 call reada(controlcard,'BOXY',boxysize,100.0d0)
223 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
224 c Cutoff range for interactions
225 call reada(controlcard,"R_CUT",r_cut,15.0d0)
226 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
227 if (me.eq.king .or. .not.out1file )
228 & write (iout,*) "DISTCHAINMAX",distchainmax
230 if(me.eq.king.or..not.out1file)
231 & write (iout,'(2a)') diagmeth(kdiag),
232 & ' routine used to diagonalize matrices.'
235 c--------------------------------------------------------------------------
236 subroutine read_REMDpar
240 implicit real*8 (a-h,o-z)
242 include 'COMMON.IOUNITS'
243 include 'COMMON.TIME1'
246 include 'COMMON.LANGEVIN'
248 include 'COMMON.LANGEVIN.lang0'
250 include 'COMMON.INTERACT'
251 include 'COMMON.NAMES'
253 include 'COMMON.REMD'
254 include 'COMMON.CONTROL'
255 include 'COMMON.SETUP'
257 character*320 controlcard
258 character*3200 controlcard1
259 integer iremd_m_total
261 if(me.eq.king.or..not.out1file)
262 & write (iout,*) "REMD setup"
264 call card_concat(controlcard)
265 call readi(controlcard,"NREP",nrep,3)
266 call readi(controlcard,"NSTEX",nstex,1000)
267 call reada(controlcard,"RETMIN",retmin,10.0d0)
268 call reada(controlcard,"RETMAX",retmax,1000.0d0)
269 mremdsync=(index(controlcard,'SYNC').gt.0)
270 call readi(controlcard,"NSYN",i_sync_step,100)
271 restart1file=(index(controlcard,'REST1FILE').gt.0)
272 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
273 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
274 if(max_cache_traj_use.gt.max_cache_traj)
275 & max_cache_traj_use=max_cache_traj
276 if(me.eq.king.or..not.out1file) then
277 cd if (traj1file) then
278 crc caching is in testing - NTWX is not ignored
279 cd write (iout,*) "NTWX value is ignored"
280 cd write (iout,*) " trajectory is stored to one file by master"
281 cd write (iout,*) " before exchange at NSTEX intervals"
283 write (iout,*) "NREP= ",nrep
284 write (iout,*) "NSTEX= ",nstex
285 write (iout,*) "SYNC= ",mremdsync
286 write (iout,*) "NSYN= ",i_sync_step
287 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
290 if (index(controlcard,'TLIST').gt.0) then
292 call card_concat(controlcard1)
293 read(controlcard1,*) (remd_t(i),i=1,nrep)
294 if(me.eq.king.or..not.out1file)
295 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
298 if (index(controlcard,'MLIST').gt.0) then
300 call card_concat(controlcard1)
301 read(controlcard1,*) (remd_m(i),i=1,nrep)
302 if(me.eq.king.or..not.out1file) then
303 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
306 iremd_m_total=iremd_m_total+remd_m(i)
308 write (iout,*) 'Total number of replicas ',iremd_m_total
311 if(me.eq.king.or..not.out1file)
312 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
315 c--------------------------------------------------------------------------
316 subroutine read_MDpar
320 implicit real*8 (a-h,o-z)
322 include 'COMMON.IOUNITS'
323 include 'COMMON.TIME1'
326 include 'COMMON.LANGEVIN'
328 include 'COMMON.LANGEVIN.lang0'
330 include 'COMMON.INTERACT'
331 include 'COMMON.NAMES'
333 include 'COMMON.SETUP'
334 include 'COMMON.CONTROL'
335 include 'COMMON.SPLITELE'
337 character*320 controlcard
339 call card_concat(controlcard)
340 call readi(controlcard,"NSTEP",n_timestep,1000000)
341 call readi(controlcard,"NTWE",ntwe,100)
342 call readi(controlcard,"NTWX",ntwx,1000)
343 call reada(controlcard,"DT",d_time,1.0d-1)
344 call reada(controlcard,"DVMAX",dvmax,2.0d1)
345 call reada(controlcard,"DAMAX",damax,1.0d1)
346 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
347 call readi(controlcard,"LANG",lang,0)
348 RESPA = index(controlcard,"RESPA") .gt. 0
349 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
350 ntime_split0=ntime_split
351 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
352 ntime_split0=ntime_split
353 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
354 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
355 rest = index(controlcard,"REST").gt.0
356 tbf = index(controlcard,"TBF").gt.0
357 usampl = index(controlcard,"USAMPL").gt.0
359 mdpdb = index(controlcard,"MDPDB").gt.0
360 call reada(controlcard,"T_BATH",t_bath,300.0d0)
361 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
362 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
363 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
364 if (count_reset_moment.eq.0) count_reset_moment=1000000000
365 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
366 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
367 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
368 if (count_reset_vel.eq.0) count_reset_vel=1000000000
369 large = index(controlcard,"LARGE").gt.0
370 print_compon = index(controlcard,"PRINT_COMPON").gt.0
371 rattle = index(controlcard,"RATTLE").gt.0
372 c if performing umbrella sampling, fragments constrained are read from the fragment file
378 if(me.eq.king.or..not.out1file) then
380 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
382 write (iout,'(a)') "The units are:"
383 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
384 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
385 & " acceleration: angstrom/(48.9 fs)**2"
386 write (iout,'(a)') "energy: kcal/mol, temperature: K"
388 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
389 write (iout,'(a60,f10.5,a)')
390 & "Initial time step of numerical integration:",d_time,
392 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
394 write (iout,'(2a,i4,a)')
395 & "A-MTS algorithm used; initial time step for fast-varying",
396 & " short-range forces split into",ntime_split," steps."
397 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
398 & r_cut," lambda",rlamb
400 write (iout,'(2a,f10.5)')
401 & "Maximum acceleration threshold to reduce the time step",
402 & "/increase split number:",damax
403 write (iout,'(2a,f10.5)')
404 & "Maximum predicted energy drift to reduce the timestep",
405 & "/increase split number:",edriftmax
406 write (iout,'(a60,f10.5)')
407 & "Maximum velocity threshold to reduce velocities:",dvmax
408 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
409 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
410 if (rattle) write (iout,'(a60)')
411 & "Rattle algorithm used to constrain the virtual bonds"
415 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
416 call reada(controlcard,"RWAT",rwat,1.4d0)
417 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
418 surfarea=index(controlcard,"SURFAREA").gt.0
419 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
420 if(me.eq.king.or..not.out1file)then
421 write (iout,'(/a,$)') "Langevin dynamics calculation"
424 & " with direct integration of Langevin equations"
425 else if (lang.eq.2) then
426 write (iout,'(a/)') " with TINKER stochasic MD integrator"
427 else if (lang.eq.3) then
428 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
429 else if (lang.eq.4) then
430 write (iout,'(a/)') " in overdamped mode"
432 write (iout,'(//a,i5)')
433 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
436 write (iout,'(a60,f10.5)') "Temperature:",t_bath
437 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
438 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
439 write (iout,'(a60,f10.5)')
440 & "Scaling factor of the friction forces:",scal_fric
441 if (surfarea) write (iout,'(2a,i10,a)')
442 & "Friction coefficients will be scaled by solvent-accessible",
443 & " surface area every",reset_fricmat," steps."
445 c Calculate friction coefficients and bounds of stochastic forces
446 eta=6*pi*cPoise*etawat
447 if(me.eq.king.or..not.out1file)
448 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
450 gamp=scal_fric*(pstok+rwat)*eta
451 stdfp=dsqrt(2*Rb*t_bath/d_time)
453 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
454 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
456 if(me.eq.king.or..not.out1file)then
457 write (iout,'(/2a/)')
458 & "Radii of site types and friction coefficients and std's of",
459 & " stochastic forces of fully exposed sites"
460 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
462 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
463 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
467 if(me.eq.king.or..not.out1file)then
468 write (iout,'(a)') "Berendsen bath calculation"
469 write (iout,'(a60,f10.5)') "Temperature:",t_bath
470 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
472 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
473 & count_reset_moment," steps"
475 & write (iout,'(a,i10,a)')
476 & "Velocities will be reset at random every",count_reset_vel,
480 if(me.eq.king.or..not.out1file)
481 & write (iout,'(a31)') "Microcanonical mode calculation"
483 if(me.eq.king.or..not.out1file)then
484 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
486 write(iout,*) "MD running with constraints."
487 write(iout,*) "Equilibration time ", eq_time, " mtus."
488 write(iout,*) "Constraining ", nfrag," fragments."
489 write(iout,*) "Length of each fragment, weight and q0:"
491 write (iout,*) "Set of restraints #",iset
493 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
494 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
496 write(iout,*) "constraints between ", npair, "fragments."
497 write(iout,*) "constraint pairs, weights and q0:"
499 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
500 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
502 write(iout,*) "angle constraints within ", nfrag_back,
503 & "backbone fragments."
504 write(iout,*) "fragment, weights:"
506 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
507 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
508 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
511 iset=mod(kolor,nset)+1
514 if(me.eq.king.or..not.out1file)
515 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
518 c------------------------------------------------------------------------------
521 C Read molecular data.
523 implicit real*8 (a-h,o-z)
529 include 'COMMON.IOUNITS'
532 include 'COMMON.INTERACT'
533 include 'COMMON.LOCAL'
534 include 'COMMON.NAMES'
535 include 'COMMON.CHAIN'
536 include 'COMMON.FFIELD'
537 include 'COMMON.SBRIDGE'
538 include 'COMMON.HEADER'
539 include 'COMMON.CONTROL'
540 include 'COMMON.DBASE'
541 include 'COMMON.THREAD'
542 include 'COMMON.CONTACTS'
543 include 'COMMON.TORCNSTR'
544 include 'COMMON.TIME1'
545 include 'COMMON.BOUNDS'
547 include 'COMMON.SETUP'
548 character*4 sequence(maxres)
550 double precision x(maxvar)
551 character*256 pdbfile
552 character*320 weightcard
553 character*80 weightcard_t,ucase
554 dimension itype_pdb(maxres)
555 common /pizda/ itype_pdb
556 logical seq_comp,fail
557 double precision energia(0:n_ene)
563 C Read weights of the subsequent energy terms.
564 call card_concat(weightcard)
565 call reada(weightcard,'WLONG',wlong,1.0D0)
566 call reada(weightcard,'WSC',wsc,wlong)
567 call reada(weightcard,'WSCP',wscp,wlong)
568 call reada(weightcard,'WELEC',welec,1.0D0)
569 call reada(weightcard,'WVDWPP',wvdwpp,welec)
570 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
571 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
572 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
573 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
574 call reada(weightcard,'WTURN3',wturn3,1.0D0)
575 call reada(weightcard,'WTURN4',wturn4,1.0D0)
576 call reada(weightcard,'WTURN6',wturn6,1.0D0)
577 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
578 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
579 call reada(weightcard,'WBOND',wbond,1.0D0)
580 call reada(weightcard,'WTOR',wtor,1.0D0)
581 call reada(weightcard,'WTORD',wtor_d,1.0D0)
582 call reada(weightcard,'WANG',wang,1.0D0)
583 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
584 call reada(weightcard,'SCAL14',scal14,0.4D0)
585 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
586 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
587 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
588 call reada(weightcard,'TEMP0',temp0,300.0d0)
589 if (index(weightcard,'SOFT').gt.0) ipot=6
590 C 12/1/95 Added weight for the multi-body term WCORR
591 call reada(weightcard,'WCORRH',wcorr,1.0D0)
592 if (wcorr4.gt.0.0d0) wcorr=wcorr4
612 if(me.eq.king.or..not.out1file)
613 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
614 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
616 10 format (/'Energy-term weights (unscaled):'//
617 & 'WSCC= ',f10.6,' (SC-SC)'/
618 & 'WSCP= ',f10.6,' (SC-p)'/
619 & 'WELEC= ',f10.6,' (p-p electr)'/
620 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
621 & 'WBOND= ',f10.6,' (stretching)'/
622 & 'WANG= ',f10.6,' (bending)'/
623 & 'WSCLOC= ',f10.6,' (SC local)'/
624 & 'WTOR= ',f10.6,' (torsional)'/
625 & 'WTORD= ',f10.6,' (double torsional)'/
626 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
627 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
628 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
629 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
630 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
631 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
632 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
633 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
634 & 'WTURN6= ',f10.6,' (turns, 6th order)')
635 if(me.eq.king.or..not.out1file)then
636 if (wcorr4.gt.0.0d0) then
637 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
638 & 'between contact pairs of peptide groups'
639 write (iout,'(2(a,f5.3/))')
640 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
641 & 'Range of quenching the correlation terms:',2*delt_corr
642 else if (wcorr.gt.0.0d0) then
643 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
644 & 'between contact pairs of peptide groups'
646 write (iout,'(a,f8.3)')
647 & 'Scaling factor of 1,4 SC-p interactions:',scal14
648 write (iout,'(a,f8.3)')
649 & 'General scaling factor of SC-p interactions:',scalscp
651 r0_corr=cutoff_corr-delt_corr
653 aad(i,1)=scalscp*aad(i,1)
654 aad(i,2)=scalscp*aad(i,2)
655 bad(i,1)=scalscp*bad(i,1)
656 bad(i,2)=scalscp*bad(i,2)
658 call rescale_weights(t_bath)
659 if(me.eq.king.or..not.out1file)
660 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
661 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
663 22 format (/'Energy-term weights (scaled):'//
664 & 'WSCC= ',f10.6,' (SC-SC)'/
665 & 'WSCP= ',f10.6,' (SC-p)'/
666 & 'WELEC= ',f10.6,' (p-p electr)'/
667 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
668 & 'WBOND= ',f10.6,' (stretching)'/
669 & 'WANG= ',f10.6,' (bending)'/
670 & 'WSCLOC= ',f10.6,' (SC local)'/
671 & 'WTOR= ',f10.6,' (torsional)'/
672 & 'WTORD= ',f10.6,' (double torsional)'/
673 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
674 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
675 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
676 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
677 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
678 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
679 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
680 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
681 & 'WTURN6= ',f10.6,' (turns, 6th order)')
682 if(me.eq.king.or..not.out1file)
683 & write (iout,*) "Reference temperature for weights calculation:",
685 call reada(weightcard,"D0CM",d0cm,3.78d0)
686 call reada(weightcard,"AKCM",akcm,15.1d0)
687 call reada(weightcard,"AKTH",akth,11.0d0)
688 call reada(weightcard,"AKCT",akct,12.0d0)
689 call reada(weightcard,"V1SS",v1ss,-1.08d0)
690 call reada(weightcard,"V2SS",v2ss,7.61d0)
691 call reada(weightcard,"V3SS",v3ss,13.7d0)
692 call reada(weightcard,"EBR",ebr,-5.50D0)
693 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
695 dyn_ss_mask(i)=.false.
699 dyn_ssbond_ij(i,j)=1.0d300
702 call reada(weightcard,"HT",Ht,0.0D0)
704 ss_depth=ebr/wsc-0.25*eps(1,1)
705 Ht=Ht/wsc-0.25*eps(1,1)
706 akcm=akcm*wstrain/wsc
707 akth=akth*wstrain/wsc
708 akct=akct*wstrain/wsc
709 v1ss=v1ss*wstrain/wsc
710 v2ss=v2ss*wstrain/wsc
711 v3ss=v3ss*wstrain/wsc
713 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
716 if(me.eq.king.or..not.out1file) then
717 write (iout,*) "Parameters of the SS-bond potential:"
718 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
720 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
721 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
722 write (iout,*)" HT",Ht
723 print *,'indpdb=',indpdb,' pdbref=',pdbref
725 if (indpdb.gt.0 .or. pdbref) then
726 read(inp,'(a)') pdbfile
727 if(me.eq.king.or..not.out1file)
728 & write (iout,'(2a)') 'PDB data will be read from file ',
729 & pdbfile(:ilen(pdbfile))
730 open(ipdbin,file=pdbfile,status='old',err=33)
732 33 write (iout,'(a)') 'Error opening PDB file.'
735 c print *,'Begin reading pdb data'
737 c print *,'Finished reading pdb data'
738 if(me.eq.king.or..not.out1file)
739 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
740 & ' nstart_sup=',nstart_sup
742 itype_pdb(i)=itype(i)
746 nct=nstart_sup+nsup-1
747 call contact(.false.,ncont_ref,icont_ref,co)
750 if(me.eq.king.or..not.out1file)
751 & write(iout,*)'Adding sidechains'
755 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
758 do while (fail.and.nsi.le.maxsi)
759 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
762 if(fail) write(iout,*)'Adding sidechain failed for res ',
763 & i,' after ',nsi,' trials'
768 if (indpdb.eq.0) then
769 C Read sequence if not taken from the pdb file.
771 c print *,'nres=',nres
772 if (iscode.gt.0) then
773 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
775 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
777 C Convert sequence to numeric code
779 itype(i)=rescode(i,sequence(i),iscode)
781 C Assign initial virtual bond lengths
787 vbld(i+nres)=dsc(iabs(itype(i)))
788 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
789 c write (iout,*) "i",i," itype",itype(i),
790 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
794 c print '(20i4)',(itype(i),i=1,nres)
797 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
799 if (itype(i).eq.ntyp1) then
803 else if (iabs(itype(i+1)).ne.20) then
805 else if (iabs(itype(i)).ne.20) then
812 if(me.eq.king.or..not.out1file)then
813 write (iout,*) "ITEL"
815 write (iout,*) i,itype(i),itel(i)
817 print *,'Call Read_Bridge.'
820 C 8/13/98 Set limits to generating the dihedral angles
825 read (inp,*) ndih_constr
826 if (ndih_constr.gt.0) then
828 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
829 if(me.eq.king.or..not.out1file)then
831 & 'There are',ndih_constr,' constraints on phi angles.'
833 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
837 phi0(i)=deg2rad*phi0(i)
838 drange(i)=deg2rad*drange(i)
840 if(me.eq.king.or..not.out1file)
841 & write (iout,*) 'FTORS',ftors
844 phibound(1,ii) = phi0(i)-drange(i)
845 phibound(2,ii) = phi0(i)+drange(i)
852 write (iout,'(a)') 'Boundaries in phi angle sampling:'
854 write (iout,'(a3,i5,2f10.1)')
855 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
861 cd print *,'NNT=',NNT,' NCT=',NCT
862 if (itype(1).eq.ntyp1) nnt=2
863 if (itype(nres).eq.ntyp1) nct=nct-1
865 if(me.eq.king.or..not.out1file)
866 & write (iout,'(a,i3)') 'nsup=',nsup
868 if (nsup.le.(nct-nnt+1)) then
869 do i=0,nct-nnt+1-nsup
870 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
876 & 'Error - sequences to be superposed do not match.'
879 do i=0,nsup-(nct-nnt+1)
880 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
882 nstart_sup=nstart_sup+i
888 & 'Error - sequences to be superposed do not match.'
891 if (nsup.eq.0) nsup=nct-nnt
892 if (nstart_sup.eq.0) nstart_sup=nnt
893 if (nstart_seq.eq.0) nstart_seq=nnt
894 if(me.eq.king.or..not.out1file)
895 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
896 & ' nstart_seq=',nstart_seq
898 c--- Zscore rms -------
899 if (nz_start.eq.0) nz_start=nnt
900 if (nz_end.eq.0 .and. nsup.gt.0) then
902 else if (nz_end.eq.0) then
905 if(me.eq.king.or..not.out1file)then
906 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
907 write (iout,*) 'IZ_SC=',iz_sc
909 c----------------------
912 if (.not.pdbref) then
913 call read_angles(inp,*38)
915 38 write (iout,'(a)') 'Error reading reference structure.'
917 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
918 stop 'Error reading reference structure'
922 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
932 call contact(.true.,ncont_ref,icont_ref,co)
934 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
936 if (constr_dist.gt.0) call read_dist_constr
937 write (iout,*) "After read_dist_constr nhpb",nhpb
939 if(me.eq.king.or..not.out1file)
940 & write (iout,*) 'Contact order:',co
942 if(me.eq.king.or..not.out1file)
943 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
946 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
948 if(me.eq.king.or..not.out1file)
949 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
950 & icont_ref(1,i),' ',
951 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
955 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
956 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
957 & modecalc.ne.10) then
958 C If input structure hasn't been supplied from the PDB file read or generate
960 if (iranconf.eq.0 .and. .not. extconf) then
961 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
962 & write (iout,'(a)') 'Initial geometry will be read in.'
964 read(inp,'(8f10.5)',end=36,err=36)
965 & ((c(l,k),l=1,3),k=1,nres),
966 & ((c(l,k+nres),l=1,3),k=nnt,nct)
967 write (iout,*) "Exit READ_CART"
968 write (iout,'(8f10.5)')
969 & ((c(l,k),l=1,3),k=1,nres),
970 & ((c(l,k+nres),l=1,3),k=nnt,nct)
971 call int_from_cart1(.true.)
972 write (iout,*) "Finish INT_TO_CART"
975 dc(j,i)=c(j,i+1)-c(j,i)
976 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
980 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
982 dc(j,i+nres)=c(j,i+nres)-c(j,i)
983 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
989 call read_angles(inp,*36)
992 36 write (iout,'(a)') 'Error reading angle file.'
994 call mpi_finalize( MPI_COMM_WORLD,IERR )
996 stop 'Error reading angle file.'
998 else if (extconf) then
999 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1000 & write (iout,'(a)') 'Extended chain initial geometry.'
1002 theta(i)=90d0*deg2rad
1005 phi(i)=180d0*deg2rad
1008 alph(i)=110d0*deg2rad
1011 omeg(i)=-120d0*deg2rad
1012 if (itype(i).le.0) omeg(i)=-omeg(i)
1015 if(me.eq.king.or..not.out1file)
1016 & write (iout,'(a)') 'Random-generated initial geometry.'
1020 if (me.eq.king .or. fg_rank.eq.0 .and. (
1021 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1025 call gen_rand_conf(itmp,*30)
1027 30 write (iout,*) 'Failed to generate random conformation',
1028 & ', itrial=',itrial
1029 write (*,*) 'Processor:',me,
1030 & ' Failed to generate random conformation',
1040 write (iout,'(a,i3,a)') 'Processor:',me,
1041 & ' error in generating random conformation.'
1042 write (*,'(a,i3,a)') 'Processor:',me,
1043 & ' error in generating random conformation.'
1046 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1051 & ' error in generating random conformation.'
1056 elseif (modecalc.eq.4) then
1057 read (inp,'(a)') intinname
1058 open (intin,file=intinname,status='old',err=333)
1059 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1060 & write (iout,'(a)') 'intinname',intinname
1061 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1063 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1065 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1067 stop 'Error opening angle file.'
1071 C Generate distance constraints, if the PDB structure is to be regularized.
1072 if (nthread.gt.0) then
1073 call read_threadbase
1076 if (me.eq.king .or. .not. out1file)
1078 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1079 write (iout,'(/a,i3,a)')
1080 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1081 write (iout,'(20i4)') (iss(i),i=1,ns)
1083 write(iout,*)"Running with dynamic disulfide-bond formation"
1085 write (iout,'(/a/)') 'Pre-formed links are:'
1091 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1092 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1098 if (ns.gt.0.and.dyn_ss) then
1102 forcon(i-nss)=forcon(i)
1109 dyn_ss_mask(iss(i))=.true.
1112 if (i2ndstr.gt.0) call secstrp2dihc
1113 c call geom_to_var(nvar,x)
1114 c call etotal(energia(0))
1115 c call enerprint(energia(0))
1116 c call briefout(0,etot)
1118 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1119 cd write (iout,'(a)') 'Variable list:'
1120 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1122 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1123 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1124 & 'Processor',myrank,': end reading molecular data.'
1128 c--------------------------------------------------------------------------
1129 logical function seq_comp(itypea,itypeb,length)
1131 integer length,itypea(length),itypeb(length)
1134 if (itypea(i).ne.itypeb(i)) then
1142 c-----------------------------------------------------------------------------
1143 subroutine read_bridge
1144 C Read information about disulfide bridges.
1145 implicit real*8 (a-h,o-z)
1146 include 'DIMENSIONS'
1150 include 'COMMON.IOUNITS'
1151 include 'COMMON.GEO'
1152 include 'COMMON.VAR'
1153 include 'COMMON.INTERACT'
1154 include 'COMMON.LOCAL'
1155 include 'COMMON.NAMES'
1156 include 'COMMON.CHAIN'
1157 include 'COMMON.FFIELD'
1158 include 'COMMON.SBRIDGE'
1159 include 'COMMON.HEADER'
1160 include 'COMMON.CONTROL'
1161 include 'COMMON.DBASE'
1162 include 'COMMON.THREAD'
1163 include 'COMMON.TIME1'
1164 include 'COMMON.SETUP'
1165 C Read bridging residues.
1166 read (inp,*) ns,(iss(i),i=1,ns)
1168 if(me.eq.king.or..not.out1file)
1169 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1170 C Check whether the specified bridging residues are cystines.
1172 if (itype(iss(i)).ne.1) then
1173 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1174 & 'Do you REALLY think that the residue ',
1175 & restyp(itype(iss(i))),i,
1176 & ' can form a disulfide bridge?!!!'
1177 write (*,'(2a,i3,a)')
1178 & 'Do you REALLY think that the residue ',
1179 & restyp(itype(iss(i))),i,
1180 & ' can form a disulfide bridge?!!!'
1182 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1187 C Read preformed bridges.
1189 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1191 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1194 C Check if the residues involved in bridges are in the specified list of
1195 C bridging residues.
1198 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1199 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1200 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1201 & ' contains residues present in other pairs.'
1202 write (*,'(a,i3,a)') 'Disulfide pair',i,
1203 & ' contains residues present in other pairs.'
1205 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1211 if (ihpb(i).eq.iss(j)) goto 10
1213 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1216 if (jhpb(i).eq.iss(j)) goto 20
1218 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1224 ihpb(i)=ihpb(i)+nres
1225 jhpb(i)=jhpb(i)+nres
1231 c----------------------------------------------------------------------------
1232 subroutine read_x(kanal,*)
1233 implicit real*8 (a-h,o-z)
1234 include 'DIMENSIONS'
1235 include 'COMMON.GEO'
1236 include 'COMMON.VAR'
1237 include 'COMMON.CHAIN'
1238 include 'COMMON.IOUNITS'
1239 include 'COMMON.CONTROL'
1240 include 'COMMON.LOCAL'
1241 include 'COMMON.INTERACT'
1242 c Read coordinates from input
1244 read(kanal,'(8f10.5)',end=10,err=10)
1245 & ((c(l,k),l=1,3),k=1,nres),
1246 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1249 c(j,2*nres)=c(j,nres)
1251 call int_from_cart1(.false.)
1254 dc(j,i)=c(j,i+1)-c(j,i)
1255 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1259 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1261 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1262 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1270 c----------------------------------------------------------------------------
1271 subroutine read_threadbase
1272 implicit real*8 (a-h,o-z)
1273 include 'DIMENSIONS'
1274 include 'COMMON.IOUNITS'
1275 include 'COMMON.GEO'
1276 include 'COMMON.VAR'
1277 include 'COMMON.INTERACT'
1278 include 'COMMON.LOCAL'
1279 include 'COMMON.NAMES'
1280 include 'COMMON.CHAIN'
1281 include 'COMMON.FFIELD'
1282 include 'COMMON.SBRIDGE'
1283 include 'COMMON.HEADER'
1284 include 'COMMON.CONTROL'
1285 include 'COMMON.DBASE'
1286 include 'COMMON.THREAD'
1287 include 'COMMON.TIME1'
1288 C Read pattern database for threading.
1289 read (icbase,*) nseq
1291 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1292 & nres_base(2,i),nres_base(3,i)
1293 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1295 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1296 c & nres_base(2,i),nres_base(3,i)
1297 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1301 if (weidis.eq.0.0D0) weidis=0.1D0
1310 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1311 write (iout,'(a,i5)') 'nexcl: ',nexcl
1312 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1315 c------------------------------------------------------------------------------
1316 subroutine setup_var
1317 implicit real*8 (a-h,o-z)
1318 include 'DIMENSIONS'
1319 include 'COMMON.IOUNITS'
1320 include 'COMMON.GEO'
1321 include 'COMMON.VAR'
1322 include 'COMMON.INTERACT'
1323 include 'COMMON.LOCAL'
1324 include 'COMMON.NAMES'
1325 include 'COMMON.CHAIN'
1326 include 'COMMON.FFIELD'
1327 include 'COMMON.SBRIDGE'
1328 include 'COMMON.HEADER'
1329 include 'COMMON.CONTROL'
1330 include 'COMMON.DBASE'
1331 include 'COMMON.THREAD'
1332 include 'COMMON.TIME1'
1333 C Set up variable list.
1339 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1341 ialph(i,1)=nvar+nside
1345 if (indphi.gt.0) then
1347 else if (indback.gt.0) then
1352 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1355 c----------------------------------------------------------------------------
1356 subroutine gen_dist_constr
1357 C Generate CA distance constraints.
1358 implicit real*8 (a-h,o-z)
1359 include 'DIMENSIONS'
1360 include 'COMMON.IOUNITS'
1361 include 'COMMON.GEO'
1362 include 'COMMON.VAR'
1363 include 'COMMON.INTERACT'
1364 include 'COMMON.LOCAL'
1365 include 'COMMON.NAMES'
1366 include 'COMMON.CHAIN'
1367 include 'COMMON.FFIELD'
1368 include 'COMMON.SBRIDGE'
1369 include 'COMMON.HEADER'
1370 include 'COMMON.CONTROL'
1371 include 'COMMON.DBASE'
1372 include 'COMMON.THREAD'
1373 include 'COMMON.TIME1'
1374 dimension itype_pdb(maxres)
1375 common /pizda/ itype_pdb
1377 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1378 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1379 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1381 do i=nstart_sup,nstart_sup+nsup-1
1382 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1383 cd & ' seq_pdb', restyp(itype_pdb(i))
1384 do j=i+2,nstart_sup+nsup-1
1386 ihpb(nhpb)=i+nstart_seq-nstart_sup
1387 jhpb(nhpb)=j+nstart_seq-nstart_sup
1389 dhpb(nhpb)=dist(i,j)
1392 cd write (iout,'(a)') 'Distance constraints:'
1397 cd if (ii.gt.nres) then
1402 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1403 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1404 cd & dhpb(i),forcon(i)
1408 c----------------------------------------------------------------------------
1410 implicit real*8 (a-h,o-z)
1411 include 'DIMENSIONS'
1412 include 'COMMON.MAP'
1413 include 'COMMON.IOUNITS'
1414 character*3 angid(4) /'THE','PHI','ALP','OME'/
1415 character*80 mapcard,ucase
1417 read (inp,'(a)') mapcard
1418 mapcard=ucase(mapcard)
1419 if (index(mapcard,'PHI').gt.0) then
1421 else if (index(mapcard,'THE').gt.0) then
1423 else if (index(mapcard,'ALP').gt.0) then
1425 else if (index(mapcard,'OME').gt.0) then
1428 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1429 stop 'Error - illegal variable spec in MAP card.'
1431 call readi (mapcard,'RES1',res1(imap),0)
1432 call readi (mapcard,'RES2',res2(imap),0)
1433 if (res1(imap).eq.0) then
1434 res1(imap)=res2(imap)
1435 else if (res2(imap).eq.0) then
1436 res2(imap)=res1(imap)
1438 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1440 & 'Error - illegal definition of variable group in MAP.'
1441 stop 'Error - illegal definition of variable group in MAP.'
1443 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1444 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1445 call readi(mapcard,'NSTEP',nstep(imap),0)
1446 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1448 & 'Illegal boundary and/or step size specification in MAP.'
1449 stop 'Illegal boundary and/or step size specification in MAP.'
1454 c----------------------------------------------------------------------------
1456 implicit real*8 (a-h,o-z)
1457 include 'DIMENSIONS'
1458 include 'COMMON.IOUNITS'
1459 include 'COMMON.GEO'
1460 include 'COMMON.CSA'
1461 include 'COMMON.BANK'
1462 include 'COMMON.CONTROL'
1464 character*620 mcmcard
1465 call card_concat(mcmcard)
1467 call readi(mcmcard,'NCONF',nconf,50)
1468 call readi(mcmcard,'NADD',nadd,0)
1469 call readi(mcmcard,'JSTART',jstart,1)
1470 call readi(mcmcard,'JEND',jend,1)
1471 call readi(mcmcard,'NSTMAX',nstmax,500000)
1472 call readi(mcmcard,'N0',n0,1)
1473 call readi(mcmcard,'N1',n1,6)
1474 call readi(mcmcard,'N2',n2,4)
1475 call readi(mcmcard,'N3',n3,0)
1476 call readi(mcmcard,'N4',n4,0)
1477 call readi(mcmcard,'N5',n5,0)
1478 call readi(mcmcard,'N6',n6,10)
1479 call readi(mcmcard,'N7',n7,0)
1480 call readi(mcmcard,'N8',n8,0)
1481 call readi(mcmcard,'N9',n9,0)
1482 call readi(mcmcard,'N14',n14,0)
1483 call readi(mcmcard,'N15',n15,0)
1484 call readi(mcmcard,'N16',n16,0)
1485 call readi(mcmcard,'N17',n17,0)
1486 call readi(mcmcard,'N18',n18,0)
1488 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1490 call readi(mcmcard,'NDIFF',ndiff,2)
1491 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1492 call readi(mcmcard,'IS1',is1,1)
1493 call readi(mcmcard,'IS2',is2,8)
1494 call readi(mcmcard,'NRAN0',nran0,4)
1495 call readi(mcmcard,'NRAN1',nran1,2)
1496 call readi(mcmcard,'IRR',irr,1)
1497 call readi(mcmcard,'NSEED',nseed,20)
1498 call readi(mcmcard,'NTOTAL',ntotal,10000)
1499 call reada(mcmcard,'CUT1',cut1,2.0d0)
1500 call reada(mcmcard,'CUT2',cut2,5.0d0)
1501 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1502 call readi(mcmcard,'ICMAX',icmax,3)
1503 call readi(mcmcard,'IRESTART',irestart,0)
1504 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1507 call reada(mcmcard,'DELE',dele,20.0d0)
1508 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1509 call readi(mcmcard,'IREF',iref,0)
1510 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1511 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1512 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1513 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1514 write (iout,*) "NCONF_IN",nconf_in
1517 c----------------------------------------------------------------------------
1518 cfmc subroutine mcmfread
1519 cfmc implicit real*8 (a-h,o-z)
1520 cfmc include 'DIMENSIONS'
1521 cfmc include 'COMMON.MCMF'
1522 cfmc include 'COMMON.IOUNITS'
1523 cfmc include 'COMMON.GEO'
1524 cfmc character*80 ucase
1525 cfmc character*620 mcmcard
1526 cfmc call card_concat(mcmcard)
1528 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1529 cfmc write(iout,*)'MAXRANT=',maxrant
1530 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1531 cfmc write(iout,*)'MAXFAM=',maxfam
1532 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1533 cfmc write(iout,*)'NNET1=',nnet1
1534 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1535 cfmc write(iout,*)'NNET2=',nnet2
1536 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1537 cfmc write(iout,*)'NNET3=',nnet3
1538 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1539 cfmc write(iout,*)'ILASTT=',ilastt
1540 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1541 cfmc write(iout,*)'MAXSTR=',maxstr
1542 cfmc maxstr_f=maxstr/maxfam
1543 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1544 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1545 cfmc write(iout,*)'NMCMF=',nmcmf
1546 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1547 cfmc write(iout,*)'IFOCUS=',ifocus
1548 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1549 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1550 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1551 cfmc write(iout,*)'INTPRT=',intprt
1552 cfmc call readi(mcmcard,'IPRT',iprt,100)
1553 cfmc write(iout,*)'IPRT=',iprt
1554 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1555 cfmc write(iout,*)'IMAXTR=',imaxtr
1556 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1557 cfmc write(iout,*)'MAXEVEN=',maxeven
1558 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1559 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1560 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1561 cfmc write(iout,*)'INIMIN=',inimin
1562 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1563 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1564 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1565 cfmc write(iout,*)'NTHREAD=',nthread
1566 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1567 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1568 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1569 cfmc write(iout,*)'MAXPERT=',maxpert
1570 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1571 cfmc write(iout,*)'IRMSD=',irmsd
1572 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1573 cfmc write(iout,*)'DENEMIN=',denemin
1574 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1575 cfmc write(iout,*)'RCUT1S=',rcut1s
1576 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1577 cfmc write(iout,*)'RCUT1E=',rcut1e
1578 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1579 cfmc write(iout,*)'RCUT2S=',rcut2s
1580 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1581 cfmc write(iout,*)'RCUT2E=',rcut2e
1582 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1583 cfmc write(iout,*)'DPERT1=',d_pert1
1584 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1585 cfmc write(iout,*)'DPERT1A=',d_pert1a
1586 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1587 cfmc write(iout,*)'DPERT2=',d_pert2
1588 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1589 cfmc write(iout,*)'DPERT2A=',d_pert2a
1590 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1591 cfmc write(iout,*)'DPERT2B=',d_pert2b
1592 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1593 cfmc write(iout,*)'DPERT2C=',d_pert2c
1594 cfmc d_pert1=deg2rad*d_pert1
1595 cfmc d_pert1a=deg2rad*d_pert1a
1596 cfmc d_pert2=deg2rad*d_pert2
1597 cfmc d_pert2a=deg2rad*d_pert2a
1598 cfmc d_pert2b=deg2rad*d_pert2b
1599 cfmc d_pert2c=deg2rad*d_pert2c
1600 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1601 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1602 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1603 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1604 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1605 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1606 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1607 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1608 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1609 cfmc write(iout,*)'RCUTINI=',rcutini
1610 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1611 cfmc write(iout,*)'GRAT=',grat
1612 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1613 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1617 c----------------------------------------------------------------------------
1619 implicit real*8 (a-h,o-z)
1620 include 'DIMENSIONS'
1621 include 'COMMON.MCM'
1622 include 'COMMON.MCE'
1623 include 'COMMON.IOUNITS'
1625 character*320 mcmcard
1626 call card_concat(mcmcard)
1627 call readi(mcmcard,'MAXACC',maxacc,100)
1628 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1629 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1630 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1631 call readi(mcmcard,'MAXREPM',maxrepm,200)
1632 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1633 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1634 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1635 call reada(mcmcard,'E_UP',e_up,5.0D0)
1636 call reada(mcmcard,'DELTE',delte,0.1D0)
1637 call readi(mcmcard,'NSWEEP',nsweep,5)
1638 call readi(mcmcard,'NSTEPH',nsteph,0)
1639 call readi(mcmcard,'NSTEPC',nstepc,0)
1640 call reada(mcmcard,'TMIN',tmin,298.0D0)
1641 call reada(mcmcard,'TMAX',tmax,298.0D0)
1642 call readi(mcmcard,'NWINDOW',nwindow,0)
1643 call readi(mcmcard,'PRINT_MC',print_mc,0)
1644 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1645 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1646 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1647 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1648 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1649 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1650 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1651 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1652 if (nwindow.gt.0) then
1653 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1655 winlen(i)=winend(i)-winstart(i)+1
1658 if (tmax.lt.tmin) tmax=tmin
1659 if (tmax.eq.tmin) then
1663 if (nstepc.gt.0 .and. nsteph.gt.0) then
1664 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1665 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1667 C Probabilities of different move types
1668 sumpro_type(0)=0.0D0
1669 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1670 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1671 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1672 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1673 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1674 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1675 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1677 print *,'i',i,' sumprotype',sumpro_type(i)
1678 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1679 print *,'i',i,' sumprotype',sumpro_type(i)
1683 c----------------------------------------------------------------------------
1684 subroutine read_minim
1685 implicit real*8 (a-h,o-z)
1686 include 'DIMENSIONS'
1687 include 'COMMON.MINIM'
1688 include 'COMMON.IOUNITS'
1690 character*320 minimcard
1691 call card_concat(minimcard)
1692 call readi(minimcard,'MAXMIN',maxmin,2000)
1693 call readi(minimcard,'MAXFUN',maxfun,5000)
1694 call readi(minimcard,'MINMIN',minmin,maxmin)
1695 call readi(minimcard,'MINFUN',minfun,maxmin)
1696 call reada(minimcard,'TOLF',tolf,1.0D-2)
1697 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1698 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1699 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1700 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1701 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1702 & 'Options in energy minimization:'
1703 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1704 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1705 & 'MinMin:',MinMin,' MinFun:',MinFun,
1706 & ' TolF:',TolF,' RTolF:',RTolF
1709 c----------------------------------------------------------------------------
1710 subroutine read_angles(kanal,*)
1711 implicit real*8 (a-h,o-z)
1712 include 'DIMENSIONS'
1713 include 'COMMON.GEO'
1714 include 'COMMON.VAR'
1715 include 'COMMON.CHAIN'
1716 include 'COMMON.IOUNITS'
1717 include 'COMMON.CONTROL'
1718 c Read angles from input
1720 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1721 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1722 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1723 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1726 c 9/7/01 avoid 180 deg valence angle
1727 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1729 theta(i)=deg2rad*theta(i)
1730 phi(i)=deg2rad*phi(i)
1731 alph(i)=deg2rad*alph(i)
1732 omeg(i)=deg2rad*omeg(i)
1737 c----------------------------------------------------------------------------
1738 subroutine reada(rekord,lancuch,wartosc,default)
1740 character*(*) rekord,lancuch
1741 double precision wartosc,default
1744 iread=index(rekord,lancuch)
1745 if (iread.eq.0) then
1749 iread=iread+ilen(lancuch)+1
1750 read (rekord(iread:),*,err=10,end=10) wartosc
1755 c----------------------------------------------------------------------------
1756 subroutine readi(rekord,lancuch,wartosc,default)
1758 character*(*) rekord,lancuch
1759 integer wartosc,default
1762 iread=index(rekord,lancuch)
1763 if (iread.eq.0) then
1767 iread=iread+ilen(lancuch)+1
1768 read (rekord(iread:),*,err=10,end=10) wartosc
1773 c----------------------------------------------------------------------------
1774 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1777 integer tablica(dim),default
1778 character*(*) rekord,lancuch
1785 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1786 if (iread.eq.0) return
1787 iread=iread+ilen(lancuch)+1
1788 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1791 c----------------------------------------------------------------------------
1792 subroutine multreada(rekord,lancuch,tablica,dim,default)
1795 double precision tablica(dim),default
1796 character*(*) rekord,lancuch
1803 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1804 if (iread.eq.0) return
1805 iread=iread+ilen(lancuch)+1
1806 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1809 c----------------------------------------------------------------------------
1810 subroutine openunits
1811 implicit real*8 (a-h,o-z)
1812 include 'DIMENSIONS'
1815 character*16 form,nodename
1818 include 'COMMON.SETUP'
1819 include 'COMMON.IOUNITS'
1821 include 'COMMON.CONTROL'
1822 integer lenpre,lenpot,ilen,lentmp
1824 character*3 out1file_text,ucase
1827 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1828 call getenv_loc("PREFIX",prefix)
1830 call getenv_loc("POT",pot)
1831 call getenv_loc("DIRTMP",tmpdir)
1832 call getenv_loc("CURDIR",curdir)
1833 call getenv_loc("OUT1FILE",out1file_text)
1834 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1835 out1file_text=ucase(out1file_text)
1836 if (out1file_text(1:1).eq."Y") then
1839 out1file=fg_rank.gt.0
1844 if (lentmp.gt.0) then
1845 write (*,'(80(1h!))')
1846 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1847 write (*,'(80(1h!))')
1848 write (*,*)"All output files will be on node /tmp directory."
1850 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1851 if (me.eq.king) then
1852 write (*,*) "The master node is ",nodename
1853 else if (fg_rank.eq.0) then
1854 write (*,*) "I am the CG slave node ",nodename
1856 write (*,*) "I am the FG slave node ",nodename
1859 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1860 lenpre = lentmp+lenpre+1
1862 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1863 C Get the names and open the input files
1864 #if defined(WINIFL) || defined(WINPGI)
1865 open(1,file=pref_orig(:ilen(pref_orig))//
1866 & '.inp',status='old',readonly,shared)
1867 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1868 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1869 C Get parameter filenames and open the parameter files.
1870 call getenv_loc('BONDPAR',bondname)
1871 open (ibond,file=bondname,status='old',readonly,shared)
1872 call getenv_loc('THETPAR',thetname)
1873 open (ithep,file=thetname,status='old',readonly,shared)
1874 call getenv_loc('ROTPAR',rotname)
1875 open (irotam,file=rotname,status='old',readonly,shared)
1876 call getenv_loc('TORPAR',torname)
1877 open (itorp,file=torname,status='old',readonly,shared)
1878 call getenv_loc('TORDPAR',tordname)
1879 open (itordp,file=tordname,status='old',readonly,shared)
1880 call getenv_loc('FOURIER',fouriername)
1881 open (ifourier,file=fouriername,status='old',readonly,shared)
1882 call getenv_loc('ELEPAR',elename)
1883 open (ielep,file=elename,status='old',readonly,shared)
1884 call getenv_loc('SIDEPAR',sidename)
1885 open (isidep,file=sidename,status='old',readonly,shared)
1886 #elif (defined CRAY) || (defined AIX)
1887 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1889 c print *,"Processor",myrank," opened file 1"
1890 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1891 c print *,"Processor",myrank," opened file 9"
1892 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1893 C Get parameter filenames and open the parameter files.
1894 call getenv_loc('BONDPAR',bondname)
1895 open (ibond,file=bondname,status='old',action='read')
1896 c print *,"Processor",myrank," opened file IBOND"
1897 call getenv_loc('THETPAR',thetname)
1898 open (ithep,file=thetname,status='old',action='read')
1899 c print *,"Processor",myrank," opened file ITHEP"
1900 call getenv_loc('ROTPAR',rotname)
1901 open (irotam,file=rotname,status='old',action='read')
1902 c print *,"Processor",myrank," opened file IROTAM"
1903 call getenv_loc('TORPAR',torname)
1904 open (itorp,file=torname,status='old',action='read')
1905 c print *,"Processor",myrank," opened file ITORP"
1906 call getenv_loc('TORDPAR',tordname)
1907 open (itordp,file=tordname,status='old',action='read')
1908 c print *,"Processor",myrank," opened file ITORDP"
1909 call getenv_loc('SCCORPAR',sccorname)
1910 open (isccor,file=sccorname,status='old',action='read')
1911 c print *,"Processor",myrank," opened file ISCCOR"
1912 call getenv_loc('FOURIER',fouriername)
1913 open (ifourier,file=fouriername,status='old',action='read')
1914 c print *,"Processor",myrank," opened file IFOURIER"
1915 call getenv_loc('ELEPAR',elename)
1916 open (ielep,file=elename,status='old',action='read')
1917 c print *,"Processor",myrank," opened file IELEP"
1918 call getenv_loc('SIDEPAR',sidename)
1919 open (isidep,file=sidename,status='old',action='read')
1920 c print *,"Processor",myrank," opened file ISIDEP"
1921 c print *,"Processor",myrank," opened parameter files"
1923 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1924 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1925 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1926 C Get parameter filenames and open the parameter files.
1927 call getenv_loc('BONDPAR',bondname)
1928 open (ibond,file=bondname,status='old')
1929 call getenv_loc('THETPAR',thetname)
1930 open (ithep,file=thetname,status='old')
1931 call getenv_loc('ROTPAR',rotname)
1932 open (irotam,file=rotname,status='old')
1933 call getenv_loc('TORPAR',torname)
1934 open (itorp,file=torname,status='old')
1935 call getenv_loc('TORDPAR',tordname)
1936 open (itordp,file=tordname,status='old')
1937 call getenv_loc('SCCORPAR',sccorname)
1938 open (isccor,file=sccorname,status='old')
1939 call getenv_loc('FOURIER',fouriername)
1940 open (ifourier,file=fouriername,status='old')
1941 call getenv_loc('ELEPAR',elename)
1942 open (ielep,file=elename,status='old')
1943 call getenv_loc('SIDEPAR',sidename)
1944 open (isidep,file=sidename,status='old')
1946 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1948 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1949 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1950 C Get parameter filenames and open the parameter files.
1951 call getenv_loc('BONDPAR',bondname)
1952 open (ibond,file=bondname,status='old',readonly)
1953 call getenv_loc('THETPAR',thetname)
1954 open (ithep,file=thetname,status='old',readonly)
1955 call getenv_loc('ROTPAR',rotname)
1956 open (irotam,file=rotname,status='old',readonly)
1957 call getenv_loc('TORPAR',torname)
1958 open (itorp,file=torname,status='old',readonly)
1959 call getenv_loc('TORDPAR',tordname)
1960 open (itordp,file=tordname,status='old',readonly)
1961 call getenv_loc('SCCORPAR',sccorname)
1962 open (isccor,file=sccorname,status='old',readonly)
1964 call getenv_loc('THETPARPDB',thetname_pdb)
1965 print *,"thetname_pdb ",thetname_pdb
1966 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1967 print *,ithep_pdb," opened"
1969 call getenv_loc('FOURIER',fouriername)
1970 open (ifourier,file=fouriername,status='old',readonly)
1971 call getenv_loc('ELEPAR',elename)
1972 open (ielep,file=elename,status='old',readonly)
1973 call getenv_loc('SIDEPAR',sidename)
1974 open (isidep,file=sidename,status='old',readonly)
1976 call getenv_loc('ROTPARPDB',rotname_pdb)
1977 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1982 C 8/9/01 In the newest version SCp interaction constants are read from a file
1983 C Use -DOLDSCP to use hard-coded constants instead.
1985 call getenv_loc('SCPPAR',scpname)
1986 #if defined(WINIFL) || defined(WINPGI)
1987 open (iscpp,file=scpname,status='old',readonly,shared)
1988 #elif (defined CRAY) || (defined AIX)
1989 open (iscpp,file=scpname,status='old',action='read')
1991 open (iscpp,file=scpname,status='old')
1993 open (iscpp,file=scpname,status='old',readonly)
1996 call getenv_loc('PATTERN',patname)
1997 #if defined(WINIFL) || defined(WINPGI)
1998 open (icbase,file=patname,status='old',readonly,shared)
1999 #elif (defined CRAY) || (defined AIX)
2000 open (icbase,file=patname,status='old',action='read')
2002 open (icbase,file=patname,status='old')
2004 open (icbase,file=patname,status='old',readonly)
2007 C Open output file only for CG processes
2008 c print *,"Processor",myrank," fg_rank",fg_rank
2009 if (fg_rank.eq.0) then
2011 if (nodes.eq.1) then
2014 npos = dlog10(dfloat(nodes-1))+1
2016 if (npos.lt.3) npos=3
2017 write (liczba,'(i1)') npos
2018 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2020 write (liczba,form) me
2021 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2022 & liczba(:ilen(liczba))
2023 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2025 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2027 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2028 & liczba(:ilen(liczba))//'.mol2'
2029 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2030 & liczba(:ilen(liczba))//'.stat'
2032 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2033 & //liczba(:ilen(liczba))//'.stat')
2034 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2037 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2038 & liczba(:ilen(liczba))//'.const'
2043 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2044 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2045 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2046 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2047 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2049 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2051 rest2name=prefix(:ilen(prefix))//'.rst'
2053 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2056 #if defined(AIX) || defined(PGI)
2057 if (me.eq.king .or. .not. out1file)
2058 & open(iout,file=outname,status='unknown')
2060 if (fg_rank.gt.0) then
2061 write (liczba,'(i3.3)') myrank/nfgtasks
2062 write (ll,'(bz,i3.3)') fg_rank
2063 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2068 open(igeom,file=intname,status='unknown',position='append')
2069 open(ipdb,file=pdbname,status='unknown')
2070 open(imol2,file=mol2name,status='unknown')
2071 open(istat,file=statname,status='unknown',position='append')
2073 c1out open(iout,file=outname,status='unknown')
2076 if (me.eq.king .or. .not.out1file)
2077 & open(iout,file=outname,status='unknown')
2079 if (fg_rank.gt.0) then
2080 write (liczba,'(i3.3)') myrank/nfgtasks
2081 write (ll,'(bz,i3.3)') fg_rank
2082 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2087 open(igeom,file=intname,status='unknown',access='append')
2088 open(ipdb,file=pdbname,status='unknown')
2089 open(imol2,file=mol2name,status='unknown')
2090 open(istat,file=statname,status='unknown',access='append')
2092 c1out open(iout,file=outname,status='unknown')
2095 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2096 csa_seed=prefix(:lenpre)//'.CSA.seed'
2097 csa_history=prefix(:lenpre)//'.CSA.history'
2098 csa_bank=prefix(:lenpre)//'.CSA.bank'
2099 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2100 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2101 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2102 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2103 csa_int=prefix(:lenpre)//'.int'
2104 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2105 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2106 csa_in=prefix(:lenpre)//'.CSA.in'
2107 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2110 write (iout,'(80(1h-))')
2111 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2112 write (iout,'(80(1h-))')
2113 write (iout,*) "Input file : ",
2114 & pref_orig(:ilen(pref_orig))//'.inp'
2115 write (iout,*) "Output file : ",
2116 & outname(:ilen(outname))
2118 write (iout,*) "Sidechain potential file : ",
2119 & sidename(:ilen(sidename))
2121 write (iout,*) "SCp potential file : ",
2122 & scpname(:ilen(scpname))
2124 write (iout,*) "Electrostatic potential file : ",
2125 & elename(:ilen(elename))
2126 write (iout,*) "Cumulant coefficient file : ",
2127 & fouriername(:ilen(fouriername))
2128 write (iout,*) "Torsional parameter file : ",
2129 & torname(:ilen(torname))
2130 write (iout,*) "Double torsional parameter file : ",
2131 & tordname(:ilen(tordname))
2132 write (iout,*) "SCCOR parameter file : ",
2133 & sccorname(:ilen(sccorname))
2134 write (iout,*) "Bond & inertia constant file : ",
2135 & bondname(:ilen(bondname))
2136 write (iout,*) "Bending parameter file : ",
2137 & thetname(:ilen(thetname))
2138 write (iout,*) "Rotamer parameter file : ",
2139 & rotname(:ilen(rotname))
2140 write (iout,*) "Thetpdb parameter file : ",
2141 & thetname_pdb(:ilen(thetname_pdb))
2142 write (iout,*) "Threading database : ",
2143 & patname(:ilen(patname))
2145 &write (iout,*)" DIRTMP : ",
2147 write (iout,'(80(1h-))')
2151 c----------------------------------------------------------------------------
2152 subroutine card_concat(card)
2153 implicit real*8 (a-h,o-z)
2154 include 'DIMENSIONS'
2155 include 'COMMON.IOUNITS'
2157 character*80 karta,ucase
2159 read (inp,'(a)') karta
2162 do while (karta(80:80).eq.'&')
2163 card=card(:ilen(card)+1)//karta(:79)
2164 read (inp,'(a)') karta
2167 card=card(:ilen(card)+1)//karta
2170 c----------------------------------------------------------------------------------
2172 implicit real*8 (a-h,o-z)
2173 include 'DIMENSIONS'
2174 include 'COMMON.CHAIN'
2175 include 'COMMON.IOUNITS'
2177 open(irest2,file=rest2name,status='unknown')
2178 read(irest2,*) totT,EK,potE,totE,t_bath
2180 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2183 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2186 read (irest2,*) iset
2191 c---------------------------------------------------------------------------------
2192 subroutine read_fragments
2193 implicit real*8 (a-h,o-z)
2194 include 'DIMENSIONS'
2198 include 'COMMON.SETUP'
2199 include 'COMMON.CHAIN'
2200 include 'COMMON.IOUNITS'
2202 include 'COMMON.CONTROL'
2203 read(inp,*) nset,nfrag,npair,nfrag_back
2204 if(me.eq.king.or..not.out1file)
2205 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2206 & " nfrag_back",nfrag_back
2208 read(inp,*) mset(iset)
2210 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2212 if(me.eq.king.or..not.out1file)
2213 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2214 & ifrag(2,i,iset), qinfrag(i,iset)
2217 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2219 if(me.eq.king.or..not.out1file)
2220 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2221 & ipair(2,i,iset), qinpair(i,iset)
2224 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2225 & wfrag_back(3,i,iset),
2226 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2227 if(me.eq.king.or..not.out1file)
2228 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2229 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2234 c-------------------------------------------------------------------------------
2235 subroutine read_dist_constr
2236 implicit real*8 (a-h,o-z)
2237 include 'DIMENSIONS'
2241 include 'COMMON.SETUP'
2242 include 'COMMON.CONTROL'
2243 include 'COMMON.CHAIN'
2244 include 'COMMON.IOUNITS'
2245 include 'COMMON.SBRIDGE'
2246 integer ifrag_(2,100),ipair_(2,100)
2247 double precision wfrag_(100),wpair_(100)
2248 character*500 controlcard
2249 c write (iout,*) "Calling read_dist_constr"
2250 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2252 call card_concat(controlcard)
2253 call readi(controlcard,"NFRAG",nfrag_,0)
2254 call readi(controlcard,"NPAIR",npair_,0)
2255 call readi(controlcard,"NDIST",ndist_,0)
2256 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2257 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2258 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2259 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2260 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2261 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2262 c write (iout,*) "IFRAG"
2264 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2266 c write (iout,*) "IPAIR"
2268 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2272 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2273 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2274 & ifrag_(2,i)=nstart_sup+nsup-1
2275 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2277 if (wfrag_(i).gt.0.0d0) then
2278 do j=ifrag_(1,i),ifrag_(2,i)-1
2279 do k=j+1,ifrag_(2,i)
2280 c write (iout,*) "j",j," k",k
2282 if (constr_dist.eq.1) then
2287 forcon(nhpb)=wfrag_(i)
2288 else if (constr_dist.eq.2) then
2289 if (ddjk.le.dist_cut) then
2294 forcon(nhpb)=wfrag_(i)
2301 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2304 if (.not.out1file .or. me.eq.king)
2305 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2306 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2308 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2309 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2316 if (wpair_(i).gt.0.0d0) then
2324 do j=ifrag_(1,ii),ifrag_(2,ii)
2325 do k=ifrag_(1,jj),ifrag_(2,jj)
2329 forcon(nhpb)=wpair_(i)
2330 dhpb(nhpb)=dist(j,k)
2332 if (.not.out1file .or. me.eq.king)
2333 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2334 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2336 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2337 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2344 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2345 if (forcon(nhpb+1).gt.0.0d0) then
2347 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2349 if (.not.out1file .or. me.eq.king)
2350 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2351 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2353 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2354 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2361 c-------------------------------------------------------------------------------
2363 subroutine flush(iu)
2368 subroutine flush(iu)
2373 c------------------------------------------------------------------------------
2374 subroutine copy_to_tmp(source)
2375 include "DIMENSIONS"
2376 include "COMMON.IOUNITS"
2377 character*(*) source
2378 character* 256 tmpfile
2382 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2383 inquire(file=tmpfile,exist=ex)
2385 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2386 & " to temporary directory..."
2387 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2388 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2392 c------------------------------------------------------------------------------
2393 subroutine move_from_tmp(source)
2394 include "DIMENSIONS"
2395 include "COMMON.IOUNITS"
2396 character*(*) source
2399 write (*,*) "Moving ",source(:ilen(source)),
2400 & " from temporary directory to working directory"
2401 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2402 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2405 c------------------------------------------------------------------------------
2406 subroutine random_init(seed)
2408 C Initialize random number generator
2410 implicit real*8 (a-h,o-z)
2411 include 'DIMENSIONS'
2414 logical OKRandom, prng_restart
2416 integer iseed_array(4)
2418 include 'COMMON.IOUNITS'
2419 include 'COMMON.TIME1'
2420 include 'COMMON.THREAD'
2421 include 'COMMON.SBRIDGE'
2422 include 'COMMON.CONTROL'
2423 include 'COMMON.MCM'
2424 include 'COMMON.MAP'
2425 include 'COMMON.HEADER'
2426 include 'COMMON.CSA'
2427 include 'COMMON.CHAIN'
2428 include 'COMMON.MUCA'
2430 include 'COMMON.FFIELD'
2431 include 'COMMON.SETUP'
2432 iseed=-dint(dabs(seed))
2433 if (iseed.eq.0) then
2434 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2435 & 'Random seed undefined. The program will stop.'
2436 write (*,'(/80(1h*)/20x,a/80(1h*))')
2437 & 'Random seed undefined. The program will stop.'
2439 call mpi_finalize(mpi_comm_world,ierr)
2441 stop 'Bad random seed.'
2444 if (fg_rank.eq.0) then
2448 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2449 OKRandom = prng_restart(me,iseed)
2452 tmp=65536.0d0**(4-i)
2453 iseed_array(i) = dint(seed/tmp)
2454 seed=seed-iseed_array(i)*tmp
2457 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2458 & (iseed_array(i),i=1,4)
2459 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2460 & (iseed_array(i),i=1,4)
2461 OKRandom = prng_restart(me,iseed_array)
2464 c r1 = prng_next(me)
2465 r1=ran_number(0.0D0,1.0D0)
2467 & write (iout,*) 'ran_num',r1
2468 if (r1.lt.0.0d0) OKRandom=.false.
2470 if (.not.OKRandom) then
2471 write (iout,*) 'PRNG IS NOT WORKING!!!'
2472 print *,'PRNG IS NOT WORKING!!!'
2475 call mpi_abort(mpi_comm_world,error_msg,ierr)
2478 write (iout,*) 'too many processors for parallel prng'
2479 write (*,*) 'too many processors for parallel prng'
2487 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)