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 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 if (modecalc.eq.8) then
32 inquire (file="fort.40",exist=file_exist)
33 if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
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.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
74 include 'COMMON.HEADER'
76 include 'COMMON.CHAIN'
79 include 'COMMON.FFIELD'
80 include 'COMMON.INTERACT'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call readi(controlcard,'SYM',symetr,1)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
200 lmuca=index(controlcard,'MUCA').gt.0
201 call readi(controlcard,'MUCADYN',mucadyn,0)
202 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203 if (lmuca .and. (me.eq.king .or. .not.out1file ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
209 iscode=index(controlcard,'ONE_LETTER')
210 indphi=index(controlcard,'PHI')
211 indback=index(controlcard,'BACK')
212 iranconf=index(controlcard,'RAND_CONF')
213 i2ndstr=index(controlcard,'USE_SEC_PRED')
214 gradout=index(controlcard,'GRADOUT').gt.0
215 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
217 if (me.eq.king .or. .not.out1file )
218 & write (iout,*) "DISTCHAINMAX",distchainmax
220 if(me.eq.king.or..not.out1file)
221 & write (iout,'(2a)') diagmeth(kdiag),
222 & ' routine used to diagonalize matrices.'
225 c--------------------------------------------------------------------------
226 subroutine read_REMDpar
230 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
233 include 'COMMON.TIME1'
236 include 'COMMON.LANGEVIN'
238 include 'COMMON.LANGEVIN.lang0'
240 include 'COMMON.INTERACT'
241 include 'COMMON.NAMES'
243 include 'COMMON.REMD'
244 include 'COMMON.CONTROL'
245 include 'COMMON.SETUP'
247 character*320 controlcard
248 character*3200 controlcard1
249 integer iremd_m_total
251 if(me.eq.king.or..not.out1file)
252 & write (iout,*) "REMD setup"
254 call card_concat(controlcard)
255 call readi(controlcard,"NREP",nrep,3)
256 call readi(controlcard,"NSTEX",nstex,1000)
257 call reada(controlcard,"RETMIN",retmin,10.0d0)
258 call reada(controlcard,"RETMAX",retmax,1000.0d0)
259 mremdsync=(index(controlcard,'SYNC').gt.0)
260 call readi(controlcard,"NSYN",i_sync_step,100)
261 restart1file=(index(controlcard,'REST1FILE').gt.0)
262 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
263 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
264 if(max_cache_traj_use.gt.max_cache_traj)
265 & max_cache_traj_use=max_cache_traj
266 if(me.eq.king.or..not.out1file) then
267 cd if (traj1file) then
268 crc caching is in testing - NTWX is not ignored
269 cd write (iout,*) "NTWX value is ignored"
270 cd write (iout,*) " trajectory is stored to one file by master"
271 cd write (iout,*) " before exchange at NSTEX intervals"
273 write (iout,*) "NREP= ",nrep
274 write (iout,*) "NSTEX= ",nstex
275 write (iout,*) "SYNC= ",mremdsync
276 write (iout,*) "NSYN= ",i_sync_step
277 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
280 if (index(controlcard,'TLIST').gt.0) then
282 call card_concat(controlcard1)
283 read(controlcard1,*) (remd_t(i),i=1,nrep)
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
288 if (index(controlcard,'MLIST').gt.0) then
290 call card_concat(controlcard1)
291 read(controlcard1,*) (remd_m(i),i=1,nrep)
292 if(me.eq.king.or..not.out1file) then
293 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
296 iremd_m_total=iremd_m_total+remd_m(i)
298 write (iout,*) 'Total number of replicas ',iremd_m_total
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
305 c--------------------------------------------------------------------------
306 subroutine read_MDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.SETUP'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SPLITELE'
327 character*320 controlcard
329 call card_concat(controlcard)
330 call readi(controlcard,"NSTEP",n_timestep,1000000)
331 call readi(controlcard,"NTWE",ntwe,100)
332 call readi(controlcard,"NTWX",ntwx,1000)
333 call reada(controlcard,"DT",d_time,1.0d-1)
334 call reada(controlcard,"DVMAX",dvmax,2.0d1)
335 call reada(controlcard,"DAMAX",damax,1.0d1)
336 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
337 call readi(controlcard,"LANG",lang,0)
338 RESPA = index(controlcard,"RESPA") .gt. 0
339 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
340 ntime_split0=ntime_split
341 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
342 ntime_split0=ntime_split
343 call reada(controlcard,"R_CUT",r_cut,2.0d0)
344 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
345 rest = index(controlcard,"REST").gt.0
346 tbf = index(controlcard,"TBF").gt.0
347 usampl = index(controlcard,"USAMPL").gt.0
349 mdpdb = index(controlcard,"MDPDB").gt.0
350 call reada(controlcard,"T_BATH",t_bath,300.0d0)
351 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
352 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
353 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
354 if (count_reset_moment.eq.0) count_reset_moment=1000000000
355 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
356 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
357 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
358 if (count_reset_vel.eq.0) count_reset_vel=1000000000
359 large = index(controlcard,"LARGE").gt.0
360 print_compon = index(controlcard,"PRINT_COMPON").gt.0
361 rattle = index(controlcard,"RATTLE").gt.0
362 c if performing umbrella sampling, fragments constrained are read from the fragment file
368 if(me.eq.king.or..not.out1file) then
370 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
372 write (iout,'(a)') "The units are:"
373 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
374 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
375 & " acceleration: angstrom/(48.9 fs)**2"
376 write (iout,'(a)') "energy: kcal/mol, temperature: K"
378 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
379 write (iout,'(a60,f10.5,a)')
380 & "Initial time step of numerical integration:",d_time,
382 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
384 write (iout,'(2a,i4,a)')
385 & "A-MTS algorithm used; initial time step for fast-varying",
386 & " short-range forces split into",ntime_split," steps."
387 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
388 & r_cut," lambda",rlamb
390 write (iout,'(2a,f10.5)')
391 & "Maximum acceleration threshold to reduce the time step",
392 & "/increase split number:",damax
393 write (iout,'(2a,f10.5)')
394 & "Maximum predicted energy drift to reduce the timestep",
395 & "/increase split number:",edriftmax
396 write (iout,'(a60,f10.5)')
397 & "Maximum velocity threshold to reduce velocities:",dvmax
398 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
399 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
400 if (rattle) write (iout,'(a60)')
401 & "Rattle algorithm used to constrain the virtual bonds"
405 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
406 call reada(controlcard,"RWAT",rwat,1.4d0)
407 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
408 surfarea=index(controlcard,"SURFAREA").gt.0
409 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
410 if(me.eq.king.or..not.out1file)then
411 write (iout,'(/a,$)') "Langevin dynamics calculation"
414 & " with direct integration of Langevin equations"
415 else if (lang.eq.2) then
416 write (iout,'(a/)') " with TINKER stochasic MD integrator"
417 else if (lang.eq.3) then
418 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
419 else if (lang.eq.4) then
420 write (iout,'(a/)') " in overdamped mode"
422 write (iout,'(//a,i5)')
423 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
426 write (iout,'(a60,f10.5)') "Temperature:",t_bath
427 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
428 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
429 write (iout,'(a60,f10.5)')
430 & "Scaling factor of the friction forces:",scal_fric
431 if (surfarea) write (iout,'(2a,i10,a)')
432 & "Friction coefficients will be scaled by solvent-accessible",
433 & " surface area every",reset_fricmat," steps."
435 c Calculate friction coefficients and bounds of stochastic forces
436 eta=6*pi*cPoise*etawat
437 if(me.eq.king.or..not.out1file)
438 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
440 gamp=scal_fric*(pstok+rwat)*eta
441 stdfp=dsqrt(2*Rb*t_bath/d_time)
443 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/2a/)')
448 & "Radii of site types and friction coefficients and std's of",
449 & " stochastic forces of fully exposed sites"
450 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
452 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
453 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(a)') "Berendsen bath calculation"
459 write (iout,'(a60,f10.5)') "Temperature:",t_bath
460 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
462 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
463 & count_reset_moment," steps"
465 & write (iout,'(a,i10,a)')
466 & "Velocities will be reset at random every",count_reset_vel,
470 if(me.eq.king.or..not.out1file)
471 & write (iout,'(a31)') "Microcanonical mode calculation"
473 if(me.eq.king.or..not.out1file)then
474 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
476 write(iout,*) "MD running with constraints."
477 write(iout,*) "Equilibration time ", eq_time, " mtus."
478 write(iout,*) "Constraining ", nfrag," fragments."
479 write(iout,*) "Length of each fragment, weight and q0:"
481 write (iout,*) "Set of restraints #",iset
483 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
484 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
486 write(iout,*) "constraints between ", npair, "fragments."
487 write(iout,*) "constraint pairs, weights and q0:"
489 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
490 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
492 write(iout,*) "angle constraints within ", nfrag_back,
493 & "backbone fragments."
494 write(iout,*) "fragment, weights:"
496 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
497 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
498 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
501 iset=mod(kolor,nset)+1
504 if(me.eq.king.or..not.out1file)
505 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
508 c------------------------------------------------------------------------------
511 C Read molecular data.
513 implicit real*8 (a-h,o-z)
519 include 'COMMON.IOUNITS'
522 include 'COMMON.INTERACT'
523 include 'COMMON.LOCAL'
524 include 'COMMON.NAMES'
525 include 'COMMON.CHAIN'
526 include 'COMMON.FFIELD'
527 include 'COMMON.SBRIDGE'
528 include 'COMMON.HEADER'
529 include 'COMMON.CONTROL'
530 include 'COMMON.DBASE'
531 include 'COMMON.THREAD'
532 include 'COMMON.CONTACTS'
533 include 'COMMON.TORCNSTR'
534 include 'COMMON.TIME1'
535 include 'COMMON.BOUNDS'
537 include 'COMMON.SETUP'
538 character*4 sequence(maxres)
540 double precision x(maxvar)
541 character*256 pdbfile
542 character*320 weightcard
543 character*80 weightcard_t,ucase
544 dimension itype_pdb(maxres)
545 common /pizda/ itype_pdb
546 logical seq_comp,fail
547 double precision energia(0:n_ene)
553 C Read weights of the subsequent energy terms.
554 call card_concat(weightcard)
555 call reada(weightcard,'WLONG',wlong,1.0D0)
556 call reada(weightcard,'WSC',wsc,wlong)
557 call reada(weightcard,'WSCP',wscp,wlong)
558 call reada(weightcard,'WELEC',welec,1.0D0)
559 call reada(weightcard,'WVDWPP',wvdwpp,welec)
560 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
561 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
562 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
563 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
564 call reada(weightcard,'WTURN3',wturn3,1.0D0)
565 call reada(weightcard,'WTURN4',wturn4,1.0D0)
566 call reada(weightcard,'WTURN6',wturn6,1.0D0)
567 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
568 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
569 call reada(weightcard,'WBOND',wbond,1.0D0)
570 call reada(weightcard,'WTOR',wtor,1.0D0)
571 call reada(weightcard,'WTORD',wtor_d,1.0D0)
572 call reada(weightcard,'WANG',wang,1.0D0)
573 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
574 call reada(weightcard,'SCAL14',scal14,0.4D0)
575 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
576 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
577 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
578 call reada(weightcard,'TEMP0',temp0,300.0d0)
579 if (index(weightcard,'SOFT').gt.0) ipot=6
580 C 12/1/95 Added weight for the multi-body term WCORR
581 call reada(weightcard,'WCORRH',wcorr,1.0D0)
582 if (wcorr4.gt.0.0d0) wcorr=wcorr4
602 if(me.eq.king.or..not.out1file)
603 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
604 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
606 10 format (/'Energy-term weights (unscaled):'//
607 & 'WSCC= ',f10.6,' (SC-SC)'/
608 & 'WSCP= ',f10.6,' (SC-p)'/
609 & 'WELEC= ',f10.6,' (p-p electr)'/
610 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
611 & 'WBOND= ',f10.6,' (stretching)'/
612 & 'WANG= ',f10.6,' (bending)'/
613 & 'WSCLOC= ',f10.6,' (SC local)'/
614 & 'WTOR= ',f10.6,' (torsional)'/
615 & 'WTORD= ',f10.6,' (double torsional)'/
616 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
617 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
618 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
619 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
620 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
621 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
622 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
623 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
624 & 'WTURN6= ',f10.6,' (turns, 6th order)')
625 if(me.eq.king.or..not.out1file)then
626 if (wcorr4.gt.0.0d0) then
627 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
628 & 'between contact pairs of peptide groups'
629 write (iout,'(2(a,f5.3/))')
630 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
631 & 'Range of quenching the correlation terms:',2*delt_corr
632 else if (wcorr.gt.0.0d0) then
633 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
634 & 'between contact pairs of peptide groups'
636 write (iout,'(a,f8.3)')
637 & 'Scaling factor of 1,4 SC-p interactions:',scal14
638 write (iout,'(a,f8.3)')
639 & 'General scaling factor of SC-p interactions:',scalscp
641 r0_corr=cutoff_corr-delt_corr
643 aad(i,1)=scalscp*aad(i,1)
644 aad(i,2)=scalscp*aad(i,2)
645 bad(i,1)=scalscp*bad(i,1)
646 bad(i,2)=scalscp*bad(i,2)
648 call rescale_weights(t_bath)
649 if(me.eq.king.or..not.out1file)
650 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
651 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
653 22 format (/'Energy-term weights (scaled):'//
654 & 'WSCC= ',f10.6,' (SC-SC)'/
655 & 'WSCP= ',f10.6,' (SC-p)'/
656 & 'WELEC= ',f10.6,' (p-p electr)'/
657 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
658 & 'WBOND= ',f10.6,' (stretching)'/
659 & 'WANG= ',f10.6,' (bending)'/
660 & 'WSCLOC= ',f10.6,' (SC local)'/
661 & 'WTOR= ',f10.6,' (torsional)'/
662 & 'WTORD= ',f10.6,' (double torsional)'/
663 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
664 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
665 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
666 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
667 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
668 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
669 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
670 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
671 & 'WTURN6= ',f10.6,' (turns, 6th order)')
672 if(me.eq.king.or..not.out1file)
673 & write (iout,*) "Reference temperature for weights calculation:",
675 call reada(weightcard,"D0CM",d0cm,3.78d0)
676 call reada(weightcard,"AKCM",akcm,15.1d0)
677 call reada(weightcard,"AKTH",akth,11.0d0)
678 call reada(weightcard,"AKCT",akct,12.0d0)
679 call reada(weightcard,"V1SS",v1ss,-1.08d0)
680 call reada(weightcard,"V2SS",v2ss,7.61d0)
681 call reada(weightcard,"V3SS",v3ss,13.7d0)
682 call reada(weightcard,"EBR",ebr,-5.50D0)
683 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
685 dyn_ss_mask(i)=.false.
689 dyn_ssbond_ij(i,j)=1.0d300
692 call reada(weightcard,"HT",Ht,0.0D0)
694 ss_depth=ebr/wsc-0.25*eps(1,1)
695 Ht=Ht/wsc-0.25*eps(1,1)
696 akcm=akcm*wstrain/wsc
697 akth=akth*wstrain/wsc
698 akct=akct*wstrain/wsc
699 v1ss=v1ss*wstrain/wsc
700 v2ss=v2ss*wstrain/wsc
701 v3ss=v3ss*wstrain/wsc
703 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
706 if(me.eq.king.or..not.out1file) then
707 write (iout,*) "Parameters of the SS-bond potential:"
708 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
710 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
711 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
712 write (iout,*)" HT",Ht
713 print *,'indpdb=',indpdb,' pdbref=',pdbref
715 if (indpdb.gt.0 .or. pdbref) then
716 read(inp,'(a)') pdbfile
717 if(me.eq.king.or..not.out1file)
718 & write (iout,'(2a)') 'PDB data will be read from file ',
719 & pdbfile(:ilen(pdbfile))
720 open(ipdbin,file=pdbfile,status='old',err=33)
722 33 write (iout,'(a)') 'Error opening PDB file.'
725 c print *,'Begin reading pdb data'
727 c print *,'Finished reading pdb data'
728 if(me.eq.king.or..not.out1file)
729 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
730 & ' nstart_sup=',nstart_sup
732 itype_pdb(i)=itype(i)
736 nct=nstart_sup+nsup-1
737 call contact(.false.,ncont_ref,icont_ref,co)
740 if(me.eq.king.or..not.out1file)
741 & write(iout,*)'Adding sidechains'
745 if (iti.ne.10 .and. itype(i).ne.21) then
748 do while (fail.and.nsi.le.maxsi)
749 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
752 if(fail) write(iout,*)'Adding sidechain failed for res ',
753 & i,' after ',nsi,' trials'
758 if (indpdb.eq.0) then
759 C Read sequence if not taken from the pdb file.
761 c print *,'nres=',nres
762 if (iscode.gt.0) then
763 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
765 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
767 C Convert sequence to numeric code
769 itype(i)=rescode(i,sequence(i),iscode)
771 C Assign initial virtual bond lengths
777 vbld(i+nres)=dsc(itype(i))
778 vbld_inv(i+nres)=dsc_inv(itype(i))
779 c write (iout,*) "i",i," itype",itype(i),
780 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
784 c print '(20i4)',(itype(i),i=1,nres)
787 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
789 if (itype(i).eq.21) then
793 else if (itype(i+1).ne.20) then
795 else if (itype(i).ne.20) then
802 if(me.eq.king.or..not.out1file)then
803 write (iout,*) "ITEL"
805 write (iout,*) i,itype(i),itel(i)
807 print *,'Call Read_Bridge.'
810 C 8/13/98 Set limits to generating the dihedral angles
815 read (inp,*) ndih_constr
816 if (ndih_constr.gt.0) then
818 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
819 if(me.eq.king.or..not.out1file)then
821 & 'There are',ndih_constr,' constraints on phi angles.'
823 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
827 phi0(i)=deg2rad*phi0(i)
828 drange(i)=deg2rad*drange(i)
830 if(me.eq.king.or..not.out1file)
831 & write (iout,*) 'FTORS',ftors
834 phibound(1,ii) = phi0(i)-drange(i)
835 phibound(2,ii) = phi0(i)+drange(i)
842 write (iout,'(a)') 'Boundaries in phi angle sampling:'
844 write (iout,'(a3,i5,2f10.1)')
845 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
851 cd print *,'NNT=',NNT,' NCT=',NCT
852 if (itype(1).eq.21) nnt=2
853 if (itype(nres).eq.21) nct=nct-1
855 if(me.eq.king.or..not.out1file)
856 & write (iout,'(a,i3)') 'nsup=',nsup
858 if (nsup.le.(nct-nnt+1)) then
859 do i=0,nct-nnt+1-nsup
860 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
866 & 'Error - sequences to be superposed do not match.'
869 do i=0,nsup-(nct-nnt+1)
870 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
872 nstart_sup=nstart_sup+i
878 & 'Error - sequences to be superposed do not match.'
881 if (nsup.eq.0) nsup=nct-nnt
882 if (nstart_sup.eq.0) nstart_sup=nnt
883 if (nstart_seq.eq.0) nstart_seq=nnt
884 if(me.eq.king.or..not.out1file)
885 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
886 & ' nstart_seq=',nstart_seq
888 c--- Zscore rms -------
889 if (nz_start.eq.0) nz_start=nnt
890 if (nz_end.eq.0 .and. nsup.gt.0) then
892 else if (nz_end.eq.0) then
895 if(me.eq.king.or..not.out1file)then
896 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
897 write (iout,*) 'IZ_SC=',iz_sc
899 c----------------------
902 if (.not.pdbref) then
903 call read_angles(inp,*38)
905 38 write (iout,'(a)') 'Error reading reference structure.'
907 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
908 stop 'Error reading reference structure'
912 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
922 call contact(.true.,ncont_ref,icont_ref,co)
924 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
926 if (constr_dist.gt.0) call read_dist_constr
927 write (iout,*) "After read_dist_constr nhpb",nhpb
929 if(me.eq.king.or..not.out1file)
930 & write (iout,*) 'Contact order:',co
932 if(me.eq.king.or..not.out1file)
933 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
936 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
938 if(me.eq.king.or..not.out1file)
939 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
940 & icont_ref(1,i),' ',
941 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
945 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
946 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
947 & modecalc.ne.10) then
948 C If input structure hasn't been supplied from the PDB file read or generate
950 if (iranconf.eq.0 .and. .not. extconf) then
951 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
952 & write (iout,'(a)') 'Initial geometry will be read in.'
954 read(inp,'(8f10.5)',end=36,err=36)
955 & ((c(l,k),l=1,3),k=1,nres),
956 & ((c(l,k+nres),l=1,3),k=nnt,nct)
957 write (iout,*) "Exit READ_CART"
958 write (iout,'(8f10.5)')
959 & ((c(l,k),l=1,3),k=1,nres),
960 & ((c(l,k+nres),l=1,3),k=nnt,nct)
961 call int_from_cart1(.true.)
962 write (iout,*) "Finish INT_TO_CART"
965 dc(j,i)=c(j,i+1)-c(j,i)
966 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
970 if (itype(i).ne.10 .and. itype(i).ne.21) then
972 dc(j,i+nres)=c(j,i+nres)-c(j,i)
973 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
979 call read_angles(inp,*36)
982 36 write (iout,'(a)') 'Error reading angle file.'
984 call mpi_finalize( MPI_COMM_WORLD,IERR )
986 stop 'Error reading angle file.'
988 else if (extconf) then
989 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
990 & write (iout,'(a)') 'Extended chain initial geometry.'
992 theta(i)=90d0*deg2rad
998 alph(i)=110d0*deg2rad
1001 omeg(i)=-120d0*deg2rad
1004 if(me.eq.king.or..not.out1file)
1005 & write (iout,'(a)') 'Random-generated initial geometry.'
1009 if (me.eq.king .or. fg_rank.eq.0 .and. (
1010 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1014 call gen_rand_conf(itmp,*30)
1016 30 write (iout,*) 'Failed to generate random conformation',
1017 & ', itrial=',itrial
1018 write (*,*) 'Processor:',me,
1019 & ' Failed to generate random conformation',
1029 write (iout,'(a,i3,a)') 'Processor:',me,
1030 & ' error in generating random conformation.'
1031 write (*,'(a,i3,a)') 'Processor:',me,
1032 & ' error in generating random conformation.'
1035 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1040 & ' error in generating random conformation.'
1045 elseif (modecalc.eq.4) then
1046 read (inp,'(a)') intinname
1047 open (intin,file=intinname,status='old',err=333)
1048 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1049 & write (iout,'(a)') 'intinname',intinname
1050 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1052 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1054 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1056 stop 'Error opening angle file.'
1060 C Generate distance constraints, if the PDB structure is to be regularized.
1061 if (nthread.gt.0) then
1062 call read_threadbase
1065 if (me.eq.king .or. .not. out1file)
1067 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1068 write (iout,'(/a,i3,a)')
1069 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1070 write (iout,'(20i4)') (iss(i),i=1,ns)
1072 write(iout,*)"Running with dynamic disulfide-bond formation"
1074 write (iout,'(/a/)') 'Pre-formed links are:'
1080 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1081 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1087 if (ns.gt.0.and.dyn_ss) then
1091 forcon(i-nss)=forcon(i)
1098 dyn_ss_mask(iss(i))=.true.
1101 if (i2ndstr.gt.0) call secstrp2dihc
1102 c call geom_to_var(nvar,x)
1103 c call etotal(energia(0))
1104 c call enerprint(energia(0))
1105 c call briefout(0,etot)
1107 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1108 cd write (iout,'(a)') 'Variable list:'
1109 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1111 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1112 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1113 & 'Processor',myrank,': end reading molecular data.'
1117 c--------------------------------------------------------------------------
1118 logical function seq_comp(itypea,itypeb,length)
1120 integer length,itypea(length),itypeb(length)
1123 if (itypea(i).ne.itypeb(i)) then
1131 c-----------------------------------------------------------------------------
1132 subroutine read_bridge
1133 C Read information about disulfide bridges.
1134 implicit real*8 (a-h,o-z)
1135 include 'DIMENSIONS'
1139 include 'COMMON.IOUNITS'
1140 include 'COMMON.GEO'
1141 include 'COMMON.VAR'
1142 include 'COMMON.INTERACT'
1143 include 'COMMON.LOCAL'
1144 include 'COMMON.NAMES'
1145 include 'COMMON.CHAIN'
1146 include 'COMMON.FFIELD'
1147 include 'COMMON.SBRIDGE'
1148 include 'COMMON.HEADER'
1149 include 'COMMON.CONTROL'
1150 include 'COMMON.DBASE'
1151 include 'COMMON.THREAD'
1152 include 'COMMON.TIME1'
1153 include 'COMMON.SETUP'
1154 C Read bridging residues.
1155 read (inp,*) ns,(iss(i),i=1,ns)
1157 if(me.eq.king.or..not.out1file)
1158 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1159 C Check whether the specified bridging residues are cystines.
1161 if (itype(iss(i)).ne.1) then
1162 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1163 & 'Do you REALLY think that the residue ',
1164 & restyp(itype(iss(i))),i,
1165 & ' can form a disulfide bridge?!!!'
1166 write (*,'(2a,i3,a)')
1167 & 'Do you REALLY think that the residue ',
1168 & restyp(itype(iss(i))),i,
1169 & ' can form a disulfide bridge?!!!'
1171 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1176 C Read preformed bridges.
1178 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1180 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1183 C Check if the residues involved in bridges are in the specified list of
1184 C bridging residues.
1187 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1188 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1189 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1190 & ' contains residues present in other pairs.'
1191 write (*,'(a,i3,a)') 'Disulfide pair',i,
1192 & ' contains residues present in other pairs.'
1194 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1200 if (ihpb(i).eq.iss(j)) goto 10
1202 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1205 if (jhpb(i).eq.iss(j)) goto 20
1207 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1213 ihpb(i)=ihpb(i)+nres
1214 jhpb(i)=jhpb(i)+nres
1220 c----------------------------------------------------------------------------
1221 subroutine read_x(kanal,*)
1222 implicit real*8 (a-h,o-z)
1223 include 'DIMENSIONS'
1224 include 'COMMON.GEO'
1225 include 'COMMON.VAR'
1226 include 'COMMON.CHAIN'
1227 include 'COMMON.IOUNITS'
1228 include 'COMMON.CONTROL'
1229 include 'COMMON.LOCAL'
1230 include 'COMMON.INTERACT'
1231 c Read coordinates from input
1233 read(kanal,'(8f10.5)',end=10,err=10)
1234 & ((c(l,k),l=1,3),k=1,nres),
1235 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1238 c(j,2*nres)=c(j,nres)
1240 call int_from_cart1(.false.)
1243 dc(j,i)=c(j,i+1)-c(j,i)
1244 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1248 if (itype(i).ne.10 .and. itype(i).ne.21) then
1250 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1251 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1259 c----------------------------------------------------------------------------
1260 subroutine read_threadbase
1261 implicit real*8 (a-h,o-z)
1262 include 'DIMENSIONS'
1263 include 'COMMON.IOUNITS'
1264 include 'COMMON.GEO'
1265 include 'COMMON.VAR'
1266 include 'COMMON.INTERACT'
1267 include 'COMMON.LOCAL'
1268 include 'COMMON.NAMES'
1269 include 'COMMON.CHAIN'
1270 include 'COMMON.FFIELD'
1271 include 'COMMON.SBRIDGE'
1272 include 'COMMON.HEADER'
1273 include 'COMMON.CONTROL'
1274 include 'COMMON.DBASE'
1275 include 'COMMON.THREAD'
1276 include 'COMMON.TIME1'
1277 C Read pattern database for threading.
1278 read (icbase,*) nseq
1280 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1281 & nres_base(2,i),nres_base(3,i)
1282 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1284 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1285 c & nres_base(2,i),nres_base(3,i)
1286 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1290 if (weidis.eq.0.0D0) weidis=0.1D0
1299 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1300 write (iout,'(a,i5)') 'nexcl: ',nexcl
1301 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1304 c------------------------------------------------------------------------------
1305 subroutine setup_var
1306 implicit real*8 (a-h,o-z)
1307 include 'DIMENSIONS'
1308 include 'COMMON.IOUNITS'
1309 include 'COMMON.GEO'
1310 include 'COMMON.VAR'
1311 include 'COMMON.INTERACT'
1312 include 'COMMON.LOCAL'
1313 include 'COMMON.NAMES'
1314 include 'COMMON.CHAIN'
1315 include 'COMMON.FFIELD'
1316 include 'COMMON.SBRIDGE'
1317 include 'COMMON.HEADER'
1318 include 'COMMON.CONTROL'
1319 include 'COMMON.DBASE'
1320 include 'COMMON.THREAD'
1321 include 'COMMON.TIME1'
1322 C Set up variable list.
1328 if (itype(i).ne.10 .and. itype(i).ne.21) then
1330 ialph(i,1)=nvar+nside
1334 if (indphi.gt.0) then
1336 else if (indback.gt.0) then
1341 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1344 c----------------------------------------------------------------------------
1345 subroutine gen_dist_constr
1346 C Generate CA distance constraints.
1347 implicit real*8 (a-h,o-z)
1348 include 'DIMENSIONS'
1349 include 'COMMON.IOUNITS'
1350 include 'COMMON.GEO'
1351 include 'COMMON.VAR'
1352 include 'COMMON.INTERACT'
1353 include 'COMMON.LOCAL'
1354 include 'COMMON.NAMES'
1355 include 'COMMON.CHAIN'
1356 include 'COMMON.FFIELD'
1357 include 'COMMON.SBRIDGE'
1358 include 'COMMON.HEADER'
1359 include 'COMMON.CONTROL'
1360 include 'COMMON.DBASE'
1361 include 'COMMON.THREAD'
1362 include 'COMMON.TIME1'
1363 dimension itype_pdb(maxres)
1364 common /pizda/ itype_pdb
1366 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1367 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1368 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1370 do i=nstart_sup,nstart_sup+nsup-1
1371 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1372 cd & ' seq_pdb', restyp(itype_pdb(i))
1373 do j=i+2,nstart_sup+nsup-1
1375 ihpb(nhpb)=i+nstart_seq-nstart_sup
1376 jhpb(nhpb)=j+nstart_seq-nstart_sup
1378 dhpb(nhpb)=dist(i,j)
1381 cd write (iout,'(a)') 'Distance constraints:'
1386 cd if (ii.gt.nres) then
1391 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1392 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1393 cd & dhpb(i),forcon(i)
1397 c----------------------------------------------------------------------------
1399 implicit real*8 (a-h,o-z)
1400 include 'DIMENSIONS'
1401 include 'COMMON.MAP'
1402 include 'COMMON.IOUNITS'
1403 character*3 angid(4) /'THE','PHI','ALP','OME'/
1404 character*80 mapcard,ucase
1406 read (inp,'(a)') mapcard
1407 mapcard=ucase(mapcard)
1408 if (index(mapcard,'PHI').gt.0) then
1410 else if (index(mapcard,'THE').gt.0) then
1412 else if (index(mapcard,'ALP').gt.0) then
1414 else if (index(mapcard,'OME').gt.0) then
1417 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1418 stop 'Error - illegal variable spec in MAP card.'
1420 call readi (mapcard,'RES1',res1(imap),0)
1421 call readi (mapcard,'RES2',res2(imap),0)
1422 if (res1(imap).eq.0) then
1423 res1(imap)=res2(imap)
1424 else if (res2(imap).eq.0) then
1425 res2(imap)=res1(imap)
1427 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1429 & 'Error - illegal definition of variable group in MAP.'
1430 stop 'Error - illegal definition of variable group in MAP.'
1432 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1433 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1434 call readi(mapcard,'NSTEP',nstep(imap),0)
1435 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1437 & 'Illegal boundary and/or step size specification in MAP.'
1438 stop 'Illegal boundary and/or step size specification in MAP.'
1443 c----------------------------------------------------------------------------
1445 implicit real*8 (a-h,o-z)
1446 include 'DIMENSIONS'
1447 include 'COMMON.IOUNITS'
1448 include 'COMMON.GEO'
1449 include 'COMMON.CSA'
1450 include 'COMMON.BANK'
1451 include 'COMMON.CONTROL'
1453 character*620 mcmcard
1454 call card_concat(mcmcard)
1456 call readi(mcmcard,'NCONF',nconf,50)
1457 call readi(mcmcard,'NADD',nadd,0)
1458 call readi(mcmcard,'JSTART',jstart,1)
1459 call readi(mcmcard,'JEND',jend,1)
1460 call readi(mcmcard,'NSTMAX',nstmax,500000)
1461 call readi(mcmcard,'N0',n0,1)
1462 call readi(mcmcard,'N1',n1,6)
1463 call readi(mcmcard,'N2',n2,4)
1464 call readi(mcmcard,'N3',n3,0)
1465 call readi(mcmcard,'N4',n4,0)
1466 call readi(mcmcard,'N5',n5,0)
1467 call readi(mcmcard,'N6',n6,10)
1468 call readi(mcmcard,'N7',n7,0)
1469 call readi(mcmcard,'N8',n8,0)
1470 call readi(mcmcard,'N9',n9,0)
1471 call readi(mcmcard,'N14',n14,0)
1472 call readi(mcmcard,'N15',n15,0)
1473 call readi(mcmcard,'N16',n16,0)
1474 call readi(mcmcard,'N17',n17,0)
1475 call readi(mcmcard,'N18',n18,0)
1477 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1479 call readi(mcmcard,'NDIFF',ndiff,2)
1480 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1481 call readi(mcmcard,'IS1',is1,1)
1482 call readi(mcmcard,'IS2',is2,8)
1483 call readi(mcmcard,'NRAN0',nran0,4)
1484 call readi(mcmcard,'NRAN1',nran1,2)
1485 call readi(mcmcard,'IRR',irr,1)
1486 call readi(mcmcard,'NSEED',nseed,20)
1487 call readi(mcmcard,'NTOTAL',ntotal,10000)
1488 call reada(mcmcard,'CUT1',cut1,2.0d0)
1489 call reada(mcmcard,'CUT2',cut2,5.0d0)
1490 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1491 call readi(mcmcard,'ICMAX',icmax,3)
1492 call readi(mcmcard,'IRESTART',irestart,0)
1493 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1496 call reada(mcmcard,'DELE',dele,20.0d0)
1497 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1498 call readi(mcmcard,'IREF',iref,0)
1499 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1500 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1501 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1502 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1503 write (iout,*) "NCONF_IN",nconf_in
1506 c----------------------------------------------------------------------------
1507 cfmc subroutine mcmfread
1508 cfmc implicit real*8 (a-h,o-z)
1509 cfmc include 'DIMENSIONS'
1510 cfmc include 'COMMON.MCMF'
1511 cfmc include 'COMMON.IOUNITS'
1512 cfmc include 'COMMON.GEO'
1513 cfmc character*80 ucase
1514 cfmc character*620 mcmcard
1515 cfmc call card_concat(mcmcard)
1517 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1518 cfmc write(iout,*)'MAXRANT=',maxrant
1519 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1520 cfmc write(iout,*)'MAXFAM=',maxfam
1521 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1522 cfmc write(iout,*)'NNET1=',nnet1
1523 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1524 cfmc write(iout,*)'NNET2=',nnet2
1525 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1526 cfmc write(iout,*)'NNET3=',nnet3
1527 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1528 cfmc write(iout,*)'ILASTT=',ilastt
1529 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1530 cfmc write(iout,*)'MAXSTR=',maxstr
1531 cfmc maxstr_f=maxstr/maxfam
1532 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1533 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1534 cfmc write(iout,*)'NMCMF=',nmcmf
1535 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1536 cfmc write(iout,*)'IFOCUS=',ifocus
1537 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1538 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1539 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1540 cfmc write(iout,*)'INTPRT=',intprt
1541 cfmc call readi(mcmcard,'IPRT',iprt,100)
1542 cfmc write(iout,*)'IPRT=',iprt
1543 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1544 cfmc write(iout,*)'IMAXTR=',imaxtr
1545 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1546 cfmc write(iout,*)'MAXEVEN=',maxeven
1547 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1548 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1549 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1550 cfmc write(iout,*)'INIMIN=',inimin
1551 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1552 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1553 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1554 cfmc write(iout,*)'NTHREAD=',nthread
1555 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1556 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1557 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1558 cfmc write(iout,*)'MAXPERT=',maxpert
1559 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1560 cfmc write(iout,*)'IRMSD=',irmsd
1561 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1562 cfmc write(iout,*)'DENEMIN=',denemin
1563 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1564 cfmc write(iout,*)'RCUT1S=',rcut1s
1565 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1566 cfmc write(iout,*)'RCUT1E=',rcut1e
1567 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1568 cfmc write(iout,*)'RCUT2S=',rcut2s
1569 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1570 cfmc write(iout,*)'RCUT2E=',rcut2e
1571 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1572 cfmc write(iout,*)'DPERT1=',d_pert1
1573 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1574 cfmc write(iout,*)'DPERT1A=',d_pert1a
1575 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1576 cfmc write(iout,*)'DPERT2=',d_pert2
1577 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1578 cfmc write(iout,*)'DPERT2A=',d_pert2a
1579 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1580 cfmc write(iout,*)'DPERT2B=',d_pert2b
1581 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1582 cfmc write(iout,*)'DPERT2C=',d_pert2c
1583 cfmc d_pert1=deg2rad*d_pert1
1584 cfmc d_pert1a=deg2rad*d_pert1a
1585 cfmc d_pert2=deg2rad*d_pert2
1586 cfmc d_pert2a=deg2rad*d_pert2a
1587 cfmc d_pert2b=deg2rad*d_pert2b
1588 cfmc d_pert2c=deg2rad*d_pert2c
1589 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1590 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1591 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1592 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1593 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1594 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1595 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1596 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1597 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1598 cfmc write(iout,*)'RCUTINI=',rcutini
1599 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1600 cfmc write(iout,*)'GRAT=',grat
1601 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1602 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1606 c----------------------------------------------------------------------------
1608 implicit real*8 (a-h,o-z)
1609 include 'DIMENSIONS'
1610 include 'COMMON.MCM'
1611 include 'COMMON.MCE'
1612 include 'COMMON.IOUNITS'
1614 character*320 mcmcard
1615 call card_concat(mcmcard)
1616 call readi(mcmcard,'MAXACC',maxacc,100)
1617 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1618 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1619 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1620 call readi(mcmcard,'MAXREPM',maxrepm,200)
1621 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1622 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1623 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1624 call reada(mcmcard,'E_UP',e_up,5.0D0)
1625 call reada(mcmcard,'DELTE',delte,0.1D0)
1626 call readi(mcmcard,'NSWEEP',nsweep,5)
1627 call readi(mcmcard,'NSTEPH',nsteph,0)
1628 call readi(mcmcard,'NSTEPC',nstepc,0)
1629 call reada(mcmcard,'TMIN',tmin,298.0D0)
1630 call reada(mcmcard,'TMAX',tmax,298.0D0)
1631 call readi(mcmcard,'NWINDOW',nwindow,0)
1632 call readi(mcmcard,'PRINT_MC',print_mc,0)
1633 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1634 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1635 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1636 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1637 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1638 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1639 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1640 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1641 if (nwindow.gt.0) then
1642 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1644 winlen(i)=winend(i)-winstart(i)+1
1647 if (tmax.lt.tmin) tmax=tmin
1648 if (tmax.eq.tmin) then
1652 if (nstepc.gt.0 .and. nsteph.gt.0) then
1653 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1654 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1656 C Probabilities of different move types
1657 sumpro_type(0)=0.0D0
1658 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1659 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1660 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1661 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1662 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1663 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1664 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1666 print *,'i',i,' sumprotype',sumpro_type(i)
1667 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1668 print *,'i',i,' sumprotype',sumpro_type(i)
1672 c----------------------------------------------------------------------------
1673 subroutine read_minim
1674 implicit real*8 (a-h,o-z)
1675 include 'DIMENSIONS'
1676 include 'COMMON.MINIM'
1677 include 'COMMON.IOUNITS'
1679 character*320 minimcard
1680 call card_concat(minimcard)
1681 call readi(minimcard,'MAXMIN',maxmin,2000)
1682 call readi(minimcard,'MAXFUN',maxfun,5000)
1683 call readi(minimcard,'MINMIN',minmin,maxmin)
1684 call readi(minimcard,'MINFUN',minfun,maxmin)
1685 call reada(minimcard,'TOLF',tolf,1.0D-2)
1686 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1687 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1688 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1689 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1690 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1691 & 'Options in energy minimization:'
1692 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1693 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1694 & 'MinMin:',MinMin,' MinFun:',MinFun,
1695 & ' TolF:',TolF,' RTolF:',RTolF
1698 c----------------------------------------------------------------------------
1699 subroutine read_angles(kanal,*)
1700 implicit real*8 (a-h,o-z)
1701 include 'DIMENSIONS'
1702 include 'COMMON.GEO'
1703 include 'COMMON.VAR'
1704 include 'COMMON.CHAIN'
1705 include 'COMMON.IOUNITS'
1706 include 'COMMON.CONTROL'
1707 c Read angles from input
1709 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1710 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1711 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1712 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1715 c 9/7/01 avoid 180 deg valence angle
1716 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1718 theta(i)=deg2rad*theta(i)
1719 phi(i)=deg2rad*phi(i)
1720 alph(i)=deg2rad*alph(i)
1721 omeg(i)=deg2rad*omeg(i)
1726 c----------------------------------------------------------------------------
1727 subroutine reada(rekord,lancuch,wartosc,default)
1729 character*(*) rekord,lancuch
1730 double precision wartosc,default
1733 iread=index(rekord,lancuch)
1734 if (iread.eq.0) then
1738 iread=iread+ilen(lancuch)+1
1739 read (rekord(iread:),*,err=10,end=10) wartosc
1744 c----------------------------------------------------------------------------
1745 subroutine readi(rekord,lancuch,wartosc,default)
1747 character*(*) rekord,lancuch
1748 integer wartosc,default
1751 iread=index(rekord,lancuch)
1752 if (iread.eq.0) then
1756 iread=iread+ilen(lancuch)+1
1757 read (rekord(iread:),*,err=10,end=10) wartosc
1762 c----------------------------------------------------------------------------
1763 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1766 integer tablica(dim),default
1767 character*(*) rekord,lancuch
1774 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1775 if (iread.eq.0) return
1776 iread=iread+ilen(lancuch)+1
1777 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1780 c----------------------------------------------------------------------------
1781 subroutine multreada(rekord,lancuch,tablica,dim,default)
1784 double precision tablica(dim),default
1785 character*(*) rekord,lancuch
1792 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1793 if (iread.eq.0) return
1794 iread=iread+ilen(lancuch)+1
1795 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1798 c----------------------------------------------------------------------------
1799 subroutine openunits
1800 implicit real*8 (a-h,o-z)
1801 include 'DIMENSIONS'
1804 character*16 form,nodename
1807 include 'COMMON.SETUP'
1808 include 'COMMON.IOUNITS'
1810 include 'COMMON.CONTROL'
1811 integer lenpre,lenpot,ilen,lentmp
1813 character*3 out1file_text,ucase
1816 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1817 call getenv_loc("PREFIX",prefix)
1819 call getenv_loc("POT",pot)
1820 call getenv_loc("DIRTMP",tmpdir)
1821 call getenv_loc("CURDIR",curdir)
1822 call getenv_loc("OUT1FILE",out1file_text)
1823 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1824 out1file_text=ucase(out1file_text)
1825 if (out1file_text(1:1).eq."Y") then
1828 out1file=fg_rank.gt.0
1833 if (lentmp.gt.0) then
1834 write (*,'(80(1h!))')
1835 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1836 write (*,'(80(1h!))')
1837 write (*,*)"All output files will be on node /tmp directory."
1839 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1840 if (me.eq.king) then
1841 write (*,*) "The master node is ",nodename
1842 else if (fg_rank.eq.0) then
1843 write (*,*) "I am the CG slave node ",nodename
1845 write (*,*) "I am the FG slave node ",nodename
1848 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1849 lenpre = lentmp+lenpre+1
1851 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1852 C Get the names and open the input files
1853 #if defined(WINIFL) || defined(WINPGI)
1854 open(1,file=pref_orig(:ilen(pref_orig))//
1855 & '.inp',status='old',readonly,shared)
1856 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1857 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1858 C Get parameter filenames and open the parameter files.
1859 call getenv_loc('BONDPAR',bondname)
1860 open (ibond,file=bondname,status='old',readonly,shared)
1861 call getenv_loc('THETPAR',thetname)
1862 open (ithep,file=thetname,status='old',readonly,shared)
1863 call getenv_loc('ROTPAR',rotname)
1864 open (irotam,file=rotname,status='old',readonly,shared)
1865 call getenv_loc('TORPAR',torname)
1866 open (itorp,file=torname,status='old',readonly,shared)
1867 call getenv_loc('TORDPAR',tordname)
1868 open (itordp,file=tordname,status='old',readonly,shared)
1869 call getenv_loc('FOURIER',fouriername)
1870 open (ifourier,file=fouriername,status='old',readonly,shared)
1871 call getenv_loc('ELEPAR',elename)
1872 open (ielep,file=elename,status='old',readonly,shared)
1873 call getenv_loc('SIDEPAR',sidename)
1874 open (isidep,file=sidename,status='old',readonly,shared)
1875 #elif (defined CRAY) || (defined AIX)
1876 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1878 c print *,"Processor",myrank," opened file 1"
1879 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1880 c print *,"Processor",myrank," opened file 9"
1881 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1882 C Get parameter filenames and open the parameter files.
1883 call getenv_loc('BONDPAR',bondname)
1884 open (ibond,file=bondname,status='old',action='read')
1885 c print *,"Processor",myrank," opened file IBOND"
1886 call getenv_loc('THETPAR',thetname)
1887 open (ithep,file=thetname,status='old',action='read')
1888 c print *,"Processor",myrank," opened file ITHEP"
1889 call getenv_loc('ROTPAR',rotname)
1890 open (irotam,file=rotname,status='old',action='read')
1891 c print *,"Processor",myrank," opened file IROTAM"
1892 call getenv_loc('TORPAR',torname)
1893 open (itorp,file=torname,status='old',action='read')
1894 c print *,"Processor",myrank," opened file ITORP"
1895 call getenv_loc('TORDPAR',tordname)
1896 open (itordp,file=tordname,status='old',action='read')
1897 c print *,"Processor",myrank," opened file ITORDP"
1898 call getenv_loc('SCCORPAR',sccorname)
1899 open (isccor,file=sccorname,status='old',action='read')
1900 c print *,"Processor",myrank," opened file ISCCOR"
1901 call getenv_loc('FOURIER',fouriername)
1902 open (ifourier,file=fouriername,status='old',action='read')
1903 c print *,"Processor",myrank," opened file IFOURIER"
1904 call getenv_loc('ELEPAR',elename)
1905 open (ielep,file=elename,status='old',action='read')
1906 c print *,"Processor",myrank," opened file IELEP"
1907 call getenv_loc('SIDEPAR',sidename)
1908 open (isidep,file=sidename,status='old',action='read')
1909 c print *,"Processor",myrank," opened file ISIDEP"
1910 c print *,"Processor",myrank," opened parameter files"
1912 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1913 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1914 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1915 C Get parameter filenames and open the parameter files.
1916 call getenv_loc('BONDPAR',bondname)
1917 open (ibond,file=bondname,status='old')
1918 call getenv_loc('THETPAR',thetname)
1919 open (ithep,file=thetname,status='old')
1920 call getenv_loc('ROTPAR',rotname)
1921 open (irotam,file=rotname,status='old')
1922 call getenv_loc('TORPAR',torname)
1923 open (itorp,file=torname,status='old')
1924 call getenv_loc('TORDPAR',tordname)
1925 open (itordp,file=tordname,status='old')
1926 call getenv_loc('SCCORPAR',sccorname)
1927 open (isccor,file=sccorname,status='old')
1928 call getenv_loc('FOURIER',fouriername)
1929 open (ifourier,file=fouriername,status='old')
1930 call getenv_loc('ELEPAR',elename)
1931 open (ielep,file=elename,status='old')
1932 call getenv_loc('SIDEPAR',sidename)
1933 open (isidep,file=sidename,status='old')
1935 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1937 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1938 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1939 C Get parameter filenames and open the parameter files.
1940 call getenv_loc('BONDPAR',bondname)
1941 open (ibond,file=bondname,status='old',readonly)
1942 call getenv_loc('THETPAR',thetname)
1943 open (ithep,file=thetname,status='old',readonly)
1944 call getenv_loc('ROTPAR',rotname)
1945 open (irotam,file=rotname,status='old',readonly)
1946 call getenv_loc('TORPAR',torname)
1947 open (itorp,file=torname,status='old',readonly)
1948 call getenv_loc('TORDPAR',tordname)
1949 open (itordp,file=tordname,status='old',readonly)
1950 call getenv_loc('SCCORPAR',sccorname)
1951 open (isccor,file=sccorname,status='old',readonly)
1952 call getenv_loc('FOURIER',fouriername)
1953 open (ifourier,file=fouriername,status='old',readonly)
1954 call getenv_loc('ELEPAR',elename)
1955 open (ielep,file=elename,status='old',readonly)
1956 call getenv_loc('SIDEPAR',sidename)
1957 open (isidep,file=sidename,status='old',readonly)
1961 C 8/9/01 In the newest version SCp interaction constants are read from a file
1962 C Use -DOLDSCP to use hard-coded constants instead.
1964 call getenv_loc('SCPPAR',scpname)
1965 #if defined(WINIFL) || defined(WINPGI)
1966 open (iscpp,file=scpname,status='old',readonly,shared)
1967 #elif (defined CRAY) || (defined AIX)
1968 open (iscpp,file=scpname,status='old',action='read')
1970 open (iscpp,file=scpname,status='old')
1972 open (iscpp,file=scpname,status='old',readonly)
1975 call getenv_loc('PATTERN',patname)
1976 #if defined(WINIFL) || defined(WINPGI)
1977 open (icbase,file=patname,status='old',readonly,shared)
1978 #elif (defined CRAY) || (defined AIX)
1979 open (icbase,file=patname,status='old',action='read')
1981 open (icbase,file=patname,status='old')
1983 open (icbase,file=patname,status='old',readonly)
1986 C Open output file only for CG processes
1987 c print *,"Processor",myrank," fg_rank",fg_rank
1988 if (fg_rank.eq.0) then
1990 if (nodes.eq.1) then
1993 npos = dlog10(dfloat(nodes-1))+1
1995 if (npos.lt.3) npos=3
1996 write (liczba,'(i1)') npos
1997 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1999 write (liczba,form) me
2000 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2001 & liczba(:ilen(liczba))
2002 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2004 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2006 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2007 & liczba(:ilen(liczba))//'.mol2'
2008 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2009 & liczba(:ilen(liczba))//'.stat'
2011 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2012 & //liczba(:ilen(liczba))//'.stat')
2013 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2016 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2017 & liczba(:ilen(liczba))//'.const'
2022 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2023 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2024 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2025 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2026 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2028 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2030 rest2name=prefix(:ilen(prefix))//'.rst'
2032 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2035 #if defined(AIX) || defined(PGI)
2036 if (me.eq.king .or. .not. out1file)
2037 & open(iout,file=outname,status='unknown')
2039 if (fg_rank.gt.0) then
2040 write (liczba,'(i3.3)') myrank/nfgtasks
2041 write (ll,'(bz,i3.3)') fg_rank
2042 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2047 open(igeom,file=intname,status='unknown',position='append')
2048 open(ipdb,file=pdbname,status='unknown')
2049 open(imol2,file=mol2name,status='unknown')
2050 open(istat,file=statname,status='unknown',position='append')
2052 c1out open(iout,file=outname,status='unknown')
2055 if (me.eq.king .or. .not.out1file)
2056 & open(iout,file=outname,status='unknown')
2058 if (fg_rank.gt.0) then
2059 write (liczba,'(i3.3)') myrank/nfgtasks
2060 write (ll,'(bz,i3.3)') fg_rank
2061 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2066 open(igeom,file=intname,status='unknown',access='append')
2067 open(ipdb,file=pdbname,status='unknown')
2068 open(imol2,file=mol2name,status='unknown')
2069 open(istat,file=statname,status='unknown',access='append')
2071 c1out open(iout,file=outname,status='unknown')
2074 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2075 csa_seed=prefix(:lenpre)//'.CSA.seed'
2076 csa_history=prefix(:lenpre)//'.CSA.history'
2077 csa_bank=prefix(:lenpre)//'.CSA.bank'
2078 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2079 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2080 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2081 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2082 csa_int=prefix(:lenpre)//'.int'
2083 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2084 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2085 csa_in=prefix(:lenpre)//'.CSA.in'
2086 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2089 write (iout,'(80(1h-))')
2090 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2091 write (iout,'(80(1h-))')
2092 write (iout,*) "Input file : ",
2093 & pref_orig(:ilen(pref_orig))//'.inp'
2094 write (iout,*) "Output file : ",
2095 & outname(:ilen(outname))
2097 write (iout,*) "Sidechain potential file : ",
2098 & sidename(:ilen(sidename))
2100 write (iout,*) "SCp potential file : ",
2101 & scpname(:ilen(scpname))
2103 write (iout,*) "Electrostatic potential file : ",
2104 & elename(:ilen(elename))
2105 write (iout,*) "Cumulant coefficient file : ",
2106 & fouriername(:ilen(fouriername))
2107 write (iout,*) "Torsional parameter file : ",
2108 & torname(:ilen(torname))
2109 write (iout,*) "Double torsional parameter file : ",
2110 & tordname(:ilen(tordname))
2111 write (iout,*) "SCCOR parameter file : ",
2112 & sccorname(:ilen(sccorname))
2113 write (iout,*) "Bond & inertia constant file : ",
2114 & bondname(:ilen(bondname))
2115 write (iout,*) "Bending parameter file : ",
2116 & thetname(:ilen(thetname))
2117 write (iout,*) "Rotamer parameter file : ",
2118 & rotname(:ilen(rotname))
2119 write (iout,*) "Threading database : ",
2120 & patname(:ilen(patname))
2122 &write (iout,*)" DIRTMP : ",
2124 write (iout,'(80(1h-))')
2128 c----------------------------------------------------------------------------
2129 subroutine card_concat(card)
2130 implicit real*8 (a-h,o-z)
2131 include 'DIMENSIONS'
2132 include 'COMMON.IOUNITS'
2134 character*80 karta,ucase
2136 read (inp,'(a)') karta
2139 do while (karta(80:80).eq.'&')
2140 card=card(:ilen(card)+1)//karta(:79)
2141 read (inp,'(a)') karta
2144 card=card(:ilen(card)+1)//karta
2147 c----------------------------------------------------------------------------------
2149 implicit real*8 (a-h,o-z)
2150 include 'DIMENSIONS'
2151 include 'COMMON.CHAIN'
2152 include 'COMMON.IOUNITS'
2154 open(irest2,file=rest2name,status='unknown')
2155 read(irest2,*) totT,EK,potE,totE,t_bath
2157 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2160 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2163 read (irest2,*) iset
2168 c---------------------------------------------------------------------------------
2169 subroutine read_fragments
2170 implicit real*8 (a-h,o-z)
2171 include 'DIMENSIONS'
2175 include 'COMMON.SETUP'
2176 include 'COMMON.CHAIN'
2177 include 'COMMON.IOUNITS'
2179 include 'COMMON.CONTROL'
2180 read(inp,*) nset,nfrag,npair,nfrag_back
2181 if(me.eq.king.or..not.out1file)
2182 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2183 & " nfrag_back",nfrag_back
2185 read(inp,*) mset(iset)
2187 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2189 if(me.eq.king.or..not.out1file)
2190 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2191 & ifrag(2,i,iset), qinfrag(i,iset)
2194 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2196 if(me.eq.king.or..not.out1file)
2197 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2198 & ipair(2,i,iset), qinpair(i,iset)
2201 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2202 & wfrag_back(3,i,iset),
2203 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2204 if(me.eq.king.or..not.out1file)
2205 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2206 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2211 c-------------------------------------------------------------------------------
2212 subroutine read_dist_constr
2213 implicit real*8 (a-h,o-z)
2214 include 'DIMENSIONS'
2218 include 'COMMON.SETUP'
2219 include 'COMMON.CONTROL'
2220 include 'COMMON.CHAIN'
2221 include 'COMMON.IOUNITS'
2222 include 'COMMON.SBRIDGE'
2223 integer ifrag_(2,100),ipair_(2,100)
2224 double precision wfrag_(100),wpair_(100)
2225 character*500 controlcard
2226 c write (iout,*) "Calling read_dist_constr"
2227 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2229 call card_concat(controlcard)
2230 call readi(controlcard,"NFRAG",nfrag_,0)
2231 call readi(controlcard,"NPAIR",npair_,0)
2232 call readi(controlcard,"NDIST",ndist_,0)
2233 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2234 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2235 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2236 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2237 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2238 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2239 c write (iout,*) "IFRAG"
2241 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2243 c write (iout,*) "IPAIR"
2245 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2249 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2250 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2251 & ifrag_(2,i)=nstart_sup+nsup-1
2252 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2254 if (wfrag_(i).gt.0.0d0) then
2255 do j=ifrag_(1,i),ifrag_(2,i)-1
2256 do k=j+1,ifrag_(2,i)
2257 c write (iout,*) "j",j," k",k
2259 if (constr_dist.eq.1) then
2264 forcon(nhpb)=wfrag_(i)
2265 else if (constr_dist.eq.2) then
2266 if (ddjk.le.dist_cut) then
2271 forcon(nhpb)=wfrag_(i)
2278 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2281 if (.not.out1file .or. me.eq.king)
2282 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2283 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2285 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2286 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2293 if (wpair_(i).gt.0.0d0) then
2301 do j=ifrag_(1,ii),ifrag_(2,ii)
2302 do k=ifrag_(1,jj),ifrag_(2,jj)
2306 forcon(nhpb)=wpair_(i)
2307 dhpb(nhpb)=dist(j,k)
2309 if (.not.out1file .or. me.eq.king)
2310 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2311 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2313 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2314 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2321 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2322 if (forcon(nhpb+1).gt.0.0d0) then
2324 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2326 if (.not.out1file .or. me.eq.king)
2327 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2328 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2330 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2331 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2338 c-------------------------------------------------------------------------------
2340 subroutine flush(iu)
2345 subroutine flush(iu)
2350 c------------------------------------------------------------------------------
2351 subroutine copy_to_tmp(source)
2352 include "DIMENSIONS"
2353 include "COMMON.IOUNITS"
2354 character*(*) source
2355 character* 256 tmpfile
2359 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2360 inquire(file=tmpfile,exist=ex)
2362 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2363 & " to temporary directory..."
2364 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2365 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2369 c------------------------------------------------------------------------------
2370 subroutine move_from_tmp(source)
2371 include "DIMENSIONS"
2372 include "COMMON.IOUNITS"
2373 character*(*) source
2376 write (*,*) "Moving ",source(:ilen(source)),
2377 & " from temporary directory to working directory"
2378 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2379 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2382 c------------------------------------------------------------------------------
2383 subroutine random_init(seed)
2385 C Initialize random number generator
2387 implicit real*8 (a-h,o-z)
2388 include 'DIMENSIONS'
2391 logical OKRandom, prng_restart
2393 integer iseed_array(4)
2395 include 'COMMON.IOUNITS'
2396 include 'COMMON.TIME1'
2397 include 'COMMON.THREAD'
2398 include 'COMMON.SBRIDGE'
2399 include 'COMMON.CONTROL'
2400 include 'COMMON.MCM'
2401 include 'COMMON.MAP'
2402 include 'COMMON.HEADER'
2403 include 'COMMON.CSA'
2404 include 'COMMON.CHAIN'
2405 include 'COMMON.MUCA'
2407 include 'COMMON.FFIELD'
2408 include 'COMMON.SETUP'
2409 iseed=-dint(dabs(seed))
2410 if (iseed.eq.0) then
2411 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2412 & 'Random seed undefined. The program will stop.'
2413 write (*,'(/80(1h*)/20x,a/80(1h*))')
2414 & 'Random seed undefined. The program will stop.'
2416 call mpi_finalize(mpi_comm_world,ierr)
2418 stop 'Bad random seed.'
2421 if (fg_rank.eq.0) then
2425 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2426 OKRandom = prng_restart(me,iseed)
2429 tmp=65536.0d0**(4-i)
2430 iseed_array(i) = dint(seed/tmp)
2431 seed=seed-iseed_array(i)*tmp
2434 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2435 & (iseed_array(i),i=1,4)
2436 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2437 & (iseed_array(i),i=1,4)
2438 OKRandom = prng_restart(me,iseed_array)
2441 c r1 = prng_next(me)
2442 r1=ran_number(0.0D0,1.0D0)
2444 & write (iout,*) 'ran_num',r1
2445 if (r1.lt.0.0d0) OKRandom=.false.
2447 if (.not.OKRandom) then
2448 write (iout,*) 'PRNG IS NOT WORKING!!!'
2449 print *,'PRNG IS NOT WORKING!!!'
2452 call mpi_abort(mpi_comm_world,error_msg,ierr)
2455 write (iout,*) 'too many processors for parallel prng'
2456 write (*,*) 'too many processors for parallel prng'
2464 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)